OSM report 2022-2023
Marissa Dyck
2024-03-25
Before you begin
Notes
A few notes about this script.
If you are running this with the 2022-2023 data make sure you download the whole (OSM_2022-2023 GitHub repository)[https://github.com/ACMElabUvic/OSM_2022-2023] from the ACMElabUvic GitHub. This will ensure you have all the files, data, and proper folder structure you will need to run this code and associated analyses.
Also make sure you open RStudio through the R project (OSM_2022-2023.Rproj) this will automatically set your working directory to the correct place (wherever you saved the repository) and ensure you don’t have to change the file paths for some of the data.
Lastly, if you are looking to adapt this code for a future year of data, you will want to ensure you have run the ACME_camera_script_9-2-2024.R or .Rmd with your data as there is much data formatting, cleaning, and restructuring that has to be done before this code will work.
If you have question please email the most recent author, currently
Marissa A. Dyck
Postdoctoral research fellow
University of Victoria
School of Environmental Studies
Email: marissadyck17@gmail.com
Install packages
If you don’t already have the following packages installed, use the code below to install them.
install.packages('tidyverse')
install.packages('ggpubr')
install.packages('corrplot')
install.packages('Hmisc')
install.packages('glmmTMB')
install.packages('MuMIn')
install.packages('TMB', type = 'source')
install.packages('rphylopic')Load libraries
Then load the packages to your library.
library(tidyverse) # data tidying, visualization, and much more; this will load all tidyverse packages, can see complete list using tidyverse_packages()
library(ggpubr) # make modificaions to plot for publication (arrange plots)
library(PerformanceAnalytics) #Used to generate a correlation plot## Loading required package: xts## Loading required package: zoo##
## Attaching package: 'zoo'## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numeric##
## ################################### WARNING ###################################
## # We noticed you have dplyr installed. The dplyr lag() function breaks how #
## # base R's lag() function is supposed to work, which breaks lag(my_xts). #
## # #
## # Calls to lag(my_xts) that you enter or source() into this session won't #
## # work correctly. #
## # #
## # All package code is unaffected because it is protected by the R namespace #
## # mechanism. #
## # #
## # Set `options(xts.warn_dplyr_breaks_lag = FALSE)` to suppress this warning. #
## # #
## # You can use stats::lag() to make sure you're not using dplyr::lag(), or you #
## # can add conflictRules('dplyr', exclude = 'lag') to your .Rprofile to stop #
## # dplyr from breaking base R's lag() function. #
## ################################### WARNING #####################################
## Attaching package: 'xts'## The following objects are masked from 'package:dplyr':
##
## first, last##
## Attaching package: 'PerformanceAnalytics'## The following object is masked from 'package:graphics':
##
## legendlibrary(Hmisc) # used to generate histograms for all variables in data frame## Loading required package: lattice## Loading required package: survival## Loading required package: Formula##
## Attaching package: 'Hmisc'## The following objects are masked from 'package:dplyr':
##
## src, summarize## The following objects are masked from 'package:base':
##
## format.pval, unitslibrary(glmmTMB) #Constructing GLMMs## Warning in checkMatrixPackageVersion(): Package version inconsistency detected.
## TMB was built with Matrix version 1.4.1
## Current Matrix version is 1.5.3
## Please re-install 'TMB' from source using install.packages('TMB', type = 'source') or ask CRAN for a binary version of 'TMB' matching CRAN's 'Matrix' packagelibrary(MuMIn) # for model selection
library(rphylopic) # add animal silhouettes to graphs## You are using rphylopic v.1.3.0. Please remember to credit PhyloPic contributors (hint: `get_attribution()`) and cite rphylopic in your work (hint: `citation("rphylopic")`).Data
Load detection data
Read in saved and cleaned detection data from the ACME_camera_script_9-2-2024.R.
# detection data
# read in saved and cleaned detection data from the ACME_camera_script_9-2-2024.R
detections <- read_csv('data/processed/OSM_2022_ind_det.csv') %>%
# change site, species and event_id to factor
mutate_if(is.character,
as.factor)## Rows: 14102 Columns: 8
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## chr (4): array, site, species, event_id
## dbl (3): month, year, timediff
## dttm (1): datetime
##
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.Data formatting
In order to get plots that have the same formatting as last years’ report we have to do a bit of data formatting. First we need to make sure we are including the same relevant species (some were ignored for last years’ report or grouped together)
Last years report had the following species
- white-tailed deer
- snowshoe hare
- black bear
- coyote
- red squirrel
- fisher
- unknown
- moose
- lynx
- spruce grouse
- red fox
- striped skunk
- ruffed grouse
- owl
- grey wolf
- domestic dog
- cougar
- raven
- other
- mule deer
And they grouped all humans except for staff as ‘Humans’. Let’s look at the species we have in this year’s data and try to format it the same way
detections %>%
# group by array and species
group_by(array, species) %>%
summarise(n = n()) %>%
# have R print everything
print(n = nrow(.))## `summarise()` has grouped output by 'array'. You can override using the
## `.groups` argument.## # A tibble: 119 × 3
## # Groups: array [4]
## array species n
## <fct> <fct> <int>
## 1 LU01 Beaver 1
## 2 LU01 Black bear 380
## 3 LU01 Cougar 7
## 4 LU01 Coyote 581
## 5 LU01 Domestic dog 6
## 6 LU01 Fisher 111
## 7 LU01 Grey jay 14
## 8 LU01 Grey wolf 21
## 9 LU01 Human 3
## 10 LU01 Lynx 55
## 11 LU01 Moose 99
## 12 LU01 Other 1
## 13 LU01 Other birds 60
## 14 LU01 Otter 2
## 15 LU01 Owl 2
## 16 LU01 Porcupine 5
## 17 LU01 Raven 6
## 18 LU01 Red fox 50
## 19 LU01 Red squirrel 879
## 20 LU01 Ruffed grouse 14
## 21 LU01 Short-tailed weasel 5
## 22 LU01 Snowshoe hare 1443
## 23 LU01 Spruce grouse 12
## 24 LU01 Staff 71
## 25 LU01 Striped skunk 39
## 26 LU01 Unknown 210
## 27 LU01 Unknown canid 48
## 28 LU01 Unknown deer 175
## 29 LU01 Unknown mustelid 13
## 30 LU01 Unknown ungulate 8
## 31 LU01 White-tailed deer 1953
## 32 LU13 ATVer 31
## 33 LU13 Black bear 275
## 34 LU13 Caribou 3
## 35 LU13 Coyote 187
## 36 LU13 Fisher 5
## 37 LU13 Grey jay 2
## 38 LU13 Grey wolf 52
## 39 LU13 Human 2
## 40 LU13 Hunter 1
## 41 LU13 Long-tailed weasel 1
## 42 LU13 Lynx 115
## 43 LU13 Marten 27
## 44 LU13 Moose 128
## 45 LU13 Other birds 12
## 46 LU13 Owl 1
## 47 LU13 Red fox 2
## 48 LU13 Red squirrel 240
## 49 LU13 Ruffed grouse 7
## 50 LU13 Short-tailed weasel 7
## 51 LU13 Snowshoe hare 573
## 52 LU13 Spruce grouse 25
## 53 LU13 Staff 82
## 54 LU13 Striped skunk 1
## 55 LU13 Unknown 86
## 56 LU13 Unknown canid 10
## 57 LU13 Unknown deer 5
## 58 LU13 Unknown mustelid 3
## 59 LU13 White-tailed deer 86
## 60 LU13 Wolverine 8
## 61 LU15 ATVer 1
## 62 LU15 Beaver 2
## 63 LU15 Black bear 220
## 64 LU15 Canada goose 3
## 65 LU15 Caribou 51
## 66 LU15 Coyote 171
## 67 LU15 Fisher 25
## 68 LU15 Grey jay 21
## 69 LU15 Grey wolf 61
## 70 LU15 Long-tailed weasel 15
## 71 LU15 Lynx 122
## 72 LU15 Marten 63
## 73 LU15 Moose 157
## 74 LU15 Other birds 59
## 75 LU15 Otter 5
## 76 LU15 Owl 1
## 77 LU15 Red fox 39
## 78 LU15 Red squirrel 643
## 79 LU15 Ruffed grouse 11
## 80 LU15 Short-tailed weasel 7
## 81 LU15 Snowmobiler 1
## 82 LU15 Snowshoe hare 611
## 83 LU15 Spruce grouse 21
## 84 LU15 Staff 78
## 85 LU15 Unknown 98
## 86 LU15 Unknown canid 7
## 87 LU15 Unknown deer 47
## 88 LU15 Unknown mustelid 16
## 89 LU15 Unknown ungulate 5
## 90 LU15 White-tailed deer 429
## 91 LU21 Black bear 544
## 92 LU21 Canada goose 1
## 93 LU21 Caribou 16
## 94 LU21 Cougar 2
## 95 LU21 Coyote 51
## 96 LU21 Fisher 46
## 97 LU21 Grey jay 13
## 98 LU21 Grey wolf 55
## 99 LU21 Long-tailed weasel 1
## 100 LU21 Lynx 72
## 101 LU21 Marten 50
## 102 LU21 Moose 233
## 103 LU21 Other 1
## 104 LU21 Other birds 44
## 105 LU21 Owl 8
## 106 LU21 Red fox 14
## 107 LU21 Red squirrel 219
## 108 LU21 Ruffed grouse 11
## 109 LU21 Short-tailed weasel 2
## 110 LU21 Snowmobiler 6
## 111 LU21 Snowshoe hare 284
## 112 LU21 Spruce grouse 19
## 113 LU21 Staff 71
## 114 LU21 Unknown 162
## 115 LU21 Unknown canid 5
## 116 LU21 Unknown deer 65
## 117 LU21 Unknown mustelid 23
## 118 LU21 Unknown ungulate 4
## 119 LU21 White-tailed deer 839# now let's create a new data frame (tibble) to work with for the OSM figure summaries specifically
# I personally would lump all the unknown together and all the birds together but for the sake of consistency with last years' figures we will remove some entries, let's create a vector of entries to drop
species_drop <- c('Staff',
'Unknown deer',
'Unknown ungulate',
'Unknown canid',
'Unknown mustelid',
'Other birds')
# now we can create the new data frame with some changes consistent w/ choices made for 2021-2022
detections <- detections %>%
# for summarizing, lets lump all the recreational humans into "Humans"
mutate(species = recode_factor(species,
"Snowmobiler" = "Human",
"ATVer" = "Human",
'Hunter' = 'Human')) %>%
# remove species we don't want to plot
filter(!species %in% species_drop)We will also want to subset the data by landscape unit (LU) and generate a new data frame for each LU to use for plotting
# we will also want to create a data frame for each LU to plot individually
# LU1
dets_LU1 <- detections %>%
filter(array == 'LU01')
# LU13
dets_LU13 <- detections %>%
filter(array == 'LU13')
# LU15
dets_LU15 <- detections %>%
filter(array == 'LU15')
# LU21
dets_LU21 <- detections %>%
filter(array == 'LU21')ANDREW HELP
Can you make the above code into a forloop which assigns each new data frame created from subsetting as dets_LUname?
Detection plots
Detection data
Now we can apply the same data formatting for each LUs’ data frame using purrr.
We want to count the number of independent detections per species per LU to use in the detection plots
# apply the same formatting to each LU data frame using purrr map
detection_data <- list(dets_LU1,
dets_LU13,
dets_LU15,
dets_LU21) %>%
purrr::map(
~.x %>%
# group by species
group_by(species) %>%
# calculate a column with unique accounts of each species
mutate(count = n_distinct(event_id)) %>%
# keep just the columns we need
select(species, count) %>%
# keep only unique (distinct) rows so we should be left with one row per species, this helps with plotting later if you don't do it ggplot will try to count and plot each row it's annoying
distinct()) %>%
# set names of list objects
purrr::set_names('Detections LU01',
'Detections LU13',
'Detections LU15',
'Detections LU21')Detection plots
Now to graph independent detections for each LU using purrr, this avoids a TON of code repetition needed to plot each one individually
We use purrr::imap() instead of
purrr::map() because imap maintains the variable names in
our list (e.g. Detections LU01, Detections LU13, etc.) which we can then
use to title each plot.
Within purrr::imap() we just paste the code we would use
for a single ggplot since all the graphical elements (except the title
which we change with the file name [.y]) are the same
# create object detection plots which uses the detection_data list (w/ all 4 LUs)
detection_plots <- detection_data %>%
# use imap instead of map as it allows us to use .y to paste the list element names as the plot titles later
purrr::imap(
~.x %>%
# now just copy and paste the ggplot code for the detection graphs
ggplot(.,
aes(x = reorder(species, count), y = count)) +
# plot as bar graph using geom_col so we don't have to provide a y aesthetic
geom_col() +
# switch the x and y axis
coord_flip() +
# add the number of detections at the end of each bar
geom_text(aes(label = count),
color = "black",
size = 3,
hjust = -0.3,
vjust = 0.2) +
# label x and y axis with informative titles
labs(x = 'Species',
y = 'Number of Independent (30 min) Detections') +
# add title to plot with LU name the .y will take the name of whatever you named each list element in the detection_data list, so make sure this name is what you want on the ggtitle
ggtitle(.y) +
# set the theme
theme_classic() +
theme(plot.title = element_text(hjust = 0.5)))
# view plots, this will print each in it's own window so you have to scroll back in the plot viewer pane to look at each one
detection_plots## $`Detections LU01`
##
## $`Detections LU13`
##
## $`Detections LU15`
##
## $`Detections LU21`
Save detection plots
Now we want to save these plots in case we need each individual one (we will combine the detection and naive occ plots into a single figure for each LU later and use those for the OSM report, but we may want these standalone plots later so let’s save them while they are here).
We can save all the plots from the purrr iteration above using
purrr::imap. imap is used instead of map because it allows
us to retain the list object names (plot names) to paste as the file
name with the .y command.
IMPORTANT if you are using this code for a future github repo, DO NOT use .tiff as the file extension. This will cause issues when trying to push any changes to the github repo as the files are too large to meet githubs requirements
# save plots only use if needed
purrr::imap(
detection_plots,
~ggsave(.x,
file = paste0("figures/",
.y,
'.jpg'), # avoid using .tiff extension in the github repo, those files are too large to push to origin
dpi = 600,
width = 11,
height = 9,
units = 'in'))## $`Detections LU01`
## [1] "figures/Detections LU01.jpg"
##
## $`Detections LU13`
## [1] "figures/Detections LU13.jpg"
##
## $`Detections LU15`
## [1] "figures/Detections LU15.jpg"
##
## $`Detections LU21`
## [1] "figures/Detections LU21.jpg"Naive occupancy
Data
We also need to alter the detection data a bit to use for naive occupancy plots.
We will use the individual LU detection data like we did before and
use purrr::map() to apply the dame data formatting to all 4
data frames.
Here we want to calculate the total number of sites in each LU, the number of sites each species was detected at in each LU and then use both those numbers to calculate naive occupancy for each species in each LU
# First we need to alter the data frame a bit for these plots, let's create a data frame for each LU (I couldn't figure out how to do this without assigning individual data frames for each UGH)
# apply the same formatting to each data frame using purrr
occupancy_data <- list(dets_LU1,
dets_LU13,
dets_LU15,
dets_LU21) %>%
purrr::map(
~.x %>%
# calculate the total number of sites for each LU
mutate(total_sites = n_distinct(site)) %>%
# group by species to calculate the number of sites each spp occurred at
group_by(species) %>%
# add columns to count the number of sites each spp occurred at and then the naive occupancy
reframe(count = n_distinct(site),
naive_occ = count/total_sites,
ind_det = n_distinct(event_id)) %>%
# keep just the columns we need
select(species, naive_occ, ind_det) %>%
# keep only unique (distinct) rows so we should be left with one row per species, this helps with plotting
distinct()) %>%
purrr::set_names('Naive Occupancy LU01',
'Naive Occupancy LU13',
'Naive Occupancy LU15',
'Naive Occupancy LU21')Occupancy plots
Now we can graph naive occupancy for each LU using purrr, and as with the detection plots this saves a massive amount of coding using purrr to run an iteration on the data files and produce four plots at once instead of copying and pasting code for each individually
# create object occupancy_plots which uses the occupancy_data list (w/ all 4 LUs)
occupancy_plots <- occupancy_data %>%
# use imap instead of map as it allows us to use .y to paste the list element names as the plot titles later
purrr::imap(
~.x %>%
# now just copy and paste the ggplot code for the occupancy graphs
ggplot(.,
aes(x = fct_reorder(species,
ind_det), # this reorders the species so they match the order of the detection plot which makes it better for viewing when the plots are arranged together in 1 figure for each LU
y = naive_occ)) +
# plot as bars using geom_col() which uses stat = 'identity', instead of geom_bar() which will count the rows in each group and plot that instead of naive occ
geom_col() +
# flip x and y axis
coord_flip() +
# add text to end of bars that provides naive occ value
geom_text(aes(label = round(naive_occ, 2)),
size = 3,
hjust = -0.3,
vjust = 0.2) +
# relabel x and y axis and title
labs(x = 'Species',
y = 'Proportion of Sites With At Least One Detection') +
# set plot title using .y (name of list object)
ggtitle(.y) +
# set. theme elements
theme_classic()+
theme(plot.title = element_text(hjust = 0.5)))
# view plots
occupancy_plots## $`Naive Occupancy LU01`
##
## $`Naive Occupancy LU13`
##
## $`Naive Occupancy LU15`
##
## $`Naive Occupancy LU21`
Save occupancy plots
As with the detection plots, we might want these individual plots
later for something so we can use purrr::imap() to save
them to the figures folder
Again avoid using the .tiff extension in github
# save plots
purrr::imap(
occupancy_plots,
~ggsave(.x,
file = paste0("figures/",
.y,
'.jpg'), # avoid using .tiff extension in the github repo, those files are too large to push to origin
dpi = 600,
width = 11,
height = 9,
units = 'in'))## $`Naive Occupancy LU01`
## [1] "figures/Naive Occupancy LU01.jpg"
##
## $`Naive Occupancy LU13`
## [1] "figures/Naive Occupancy LU13.jpg"
##
## $`Naive Occupancy LU15`
## [1] "figures/Naive Occupancy LU15.jpg"
##
## $`Naive Occupancy LU21`
## [1] "figures/Naive Occupancy LU21.jpg"Final combined plots for report
The previous year’s report had a figure for each LU with the
detections plot on the top and the occupancy plot on the bottom so we
will recreate these for this year using ggarrange().
Unfortunately I could not figure out how to do this in purrr to reduce coding but luckily it isn’t too much repitition
# not sure I know how to do the following section in purrr just yet, but we've saved a ton of coding so far and it doesn't take much to arrange each of these individually
# LU1
# arrange the plots so each LU has a figure with detections on top and naive occ on bottom
LU1_det_occ_plots <- ggarrange(detection_plots$`Detections LU01`, occupancy_plots$`Naive Occupancy LU01`,
labels = c("A", "B"),
nrow = 2)
# view plot
LU1_det_occ_plots
# LU13
# arrange the plots so each LU has a figure with detections on top and naive occ on bottom
LU13_det_occ_plots <- ggarrange(detection_plots$`Detections LU13`, occupancy_plots$`Naive Occupancy LU13`,
labels = c("A", "B"),
nrow = 2)
# view plot
LU13_det_occ_plots
# LU15
# arrange the plots so each LU has a figure with detections on top and naive occ on bottom
LU15_det_occ_plots <- ggarrange(detection_plots$`Detections LU15`, occupancy_plots$`Naive Occupancy LU15`,
labels = c("A", "B"),
nrow = 2)
# view plot
LU15_det_occ_plots
# LU21
# arrange the plots so each LU has a figure with detections on top and naive occ on bottom
LU21_det_occ_plots <- ggarrange(detection_plots$`Detections LU21`, occupancy_plots$`Naive Occupancy LU21`,
labels = c("A", "B"),
nrow = 2)
# view plot
LU21_det_occ_plots
We can however, save all the figures again using purrr
# save all figures at once using purrr
final_det_occ_plots <- list(LU1_det_occ_plots,
LU13_det_occ_plots,
LU15_det_occ_plots,
LU21_det_occ_plots) %>%
purrr::set_names('LU01_det_occ_plots',
'LU13_det_occ_plots',
'LU15_det_occ_plots',
'LU21_det_occ_plots') %>%
purrr::imap(
~ggsave(.x,
file = paste0("figures/",
.y,
'.jpg'), # avoid using .tiff extension in the github repo, those files are too large to push to origin
dpi = 600,
width = 12,
height = 15,
units = 'in'))Analysis prep
Read in data
We need the proportional binomial data and the covariate data (from the ACME_camera_script_9-2-2024.R or .Rmd), let’s read those in now and check the structure of each
# response metric (proportional detections from the from the ACME_camera_script_9-2-2024.R or .Rmd)
prop_detections <- read_csv('data/processed/OSM_2022_proportional_detections.csv')## Rows: 152 Columns: 23
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## chr (1): site
## dbl (22): black_bear, coyote, fisher, moose, white-tailed_deer, cougar, grey...
##
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.# check variable structure
str(prop_detections)## spc_tbl_ [152 × 23] (S3: spec_tbl_df/tbl_df/tbl/data.frame)
## $ site : chr [1:152] "LU01_06" "LU01_10" "LU01_11" "LU01_13" ...
## $ black_bear : num [1:152] 7 3 4 7 8 9 4 5 7 7 ...
## $ coyote : num [1:152] 4 4 8 10 11 9 11 0 9 4 ...
## $ fisher : num [1:152] 5 3 3 3 2 1 1 1 0 3 ...
## $ moose : num [1:152] 3 2 5 9 1 0 2 4 1 0 ...
## $ white-tailed_deer : num [1:152] 12 5 12 12 13 14 15 9 12 10 ...
## $ cougar : num [1:152] 0 0 1 0 1 0 0 0 0 0 ...
## $ grey_wolf : num [1:152] 0 0 2 0 0 0 1 0 0 0 ...
## $ lynx : num [1:152] 0 0 1 0 1 1 0 0 0 2 ...
## $ red_fox : num [1:152] 0 0 2 0 0 0 0 0 4 0 ...
## $ wolverine : num [1:152] 0 0 0 0 0 0 0 0 0 0 ...
## $ caribou : num [1:152] 0 0 0 0 0 0 0 0 0 0 ...
## $ absent_black_bear : num [1:152] 5 3 8 5 4 3 8 7 5 5 ...
## $ absent_coyote : num [1:152] 10 1 6 5 3 5 4 15 6 11 ...
## $ absent_fisher : num [1:152] 9 2 11 12 12 13 14 14 15 12 ...
## $ absent_moose : num [1:152] 11 3 9 6 13 14 13 11 14 15 ...
## $ absent_white-tailed_deer: num [1:152] 2 0 2 3 1 0 0 6 3 5 ...
## $ absent_cougar : num [1:152] 14 5 13 15 13 14 15 15 15 15 ...
## $ absent_grey_wolf : num [1:152] 14 5 12 15 14 14 14 15 15 15 ...
## $ absent_lynx : num [1:152] 14 5 13 15 13 13 15 15 15 13 ...
## $ absent_red_fox : num [1:152] 14 5 12 15 14 14 15 15 11 15 ...
## $ absent_wolverine : num [1:152] 14 5 14 15 14 14 15 15 15 15 ...
## $ absent_caribou : num [1:152] 14 5 14 15 14 14 15 15 15 15 ...
## - attr(*, "spec")=
## .. cols(
## .. site = col_character(),
## .. black_bear = col_double(),
## .. coyote = col_double(),
## .. fisher = col_double(),
## .. moose = col_double(),
## .. `white-tailed_deer` = col_double(),
## .. cougar = col_double(),
## .. grey_wolf = col_double(),
## .. lynx = col_double(),
## .. red_fox = col_double(),
## .. wolverine = col_double(),
## .. caribou = col_double(),
## .. absent_black_bear = col_double(),
## .. absent_coyote = col_double(),
## .. absent_fisher = col_double(),
## .. absent_moose = col_double(),
## .. `absent_white-tailed_deer` = col_double(),
## .. absent_cougar = col_double(),
## .. absent_grey_wolf = col_double(),
## .. absent_lynx = col_double(),
## .. absent_red_fox = col_double(),
## .. absent_wolverine = col_double(),
## .. absent_caribou = col_double()
## .. )
## - attr(*, "problems")=<externalptr># model covariates (merged HFI and VEG data from the ACME_camera_script_9-2-2024.R or .Rmd)
covariates <- read_csv('data/processed/OSM_2022_covariates.csv',
# set the column types to read in correctly
col_types = cols(array = col_factor(),
camera = col_factor(),
site = col_factor(),
buff_dist = col_factor(),
.default = col_number()))
# check variable structure
str(covariates)## spc_tbl_ [3,100 × 76] (S3: spec_tbl_df/tbl_df/tbl/data.frame)
## $ array : Factor w/ 4 levels "LU13","LU15",..: 1 1 1 1 1 1 1 1 1 1 ...
## $ camera : Factor w/ 96 levels "18","15","03",..: 1 2 3 4 5 6 7 8 9 10 ...
## $ site : Factor w/ 155 levels "LU13_18","LU13_15",..: 1 2 3 4 5 6 7 8 9 10 ...
## $ buff_dist : Factor w/ 20 levels "250","500","750",..: 1 1 1 1 1 1 1 1 1 1 ...
## $ harvest_area : num [1:3100] 0 0 NA 0 0 0 0 0 0 0 ...
## $ crop : num [1:3100] 0 0 NA 0 0 0 0 0 0 0 ...
## $ well_aband : num [1:3100] 0 0 NA 0 0 0 0 0 0 0 ...
## $ well_oil : num [1:3100] 0 0 NA 0 0 0 0 0 0 0 ...
## $ trail : num [1:3100] 0 0 NA 0.5 0 ...
## $ harvest_area_white_zone : num [1:3100] 0 0 NA 0 0 0 0 0 0 0 ...
## $ conventional_seismic : num [1:3100] 0.5 0.5 NA 0.5 1 ...
## $ pipeline : num [1:3100] 0 0.5 NA 0 0 0 0.5 0 0 0 ...
## $ tame_pasture : num [1:3100] 0 0 NA 0 0 0 0 0 0 0 ...
## $ rough_pasture : num [1:3100] 0 0 NA 0 0 0 0 0 0 0 ...
## $ rural_residence : num [1:3100] 0 0 NA 0 0 0 0 0 0 0 ...
## $ transmission_line : num [1:3100] 0 0 NA 0 0 0 0 0 0 0 ...
## $ well_gas : num [1:3100] 0 0 NA 0 0 0 0 0 0 0 ...
## $ misc_oil_gas_facility : num [1:3100] 0 0 NA 0 0 0 0 0 0 0 ...
## $ clearing_unknown : num [1:3100] 0 0 NA 0 0 0 0 0 0 0 ...
## $ vegetated_edge_roads : num [1:3100] 0 0 NA 0 0 0 0 0 0 0 ...
## $ road_unimproved : num [1:3100] 0 0 NA 0 0 0 0 0 0 0 ...
## $ road_gravel_1l : num [1:3100] 0 0 NA 0 0 0 0 0 0 0 ...
## $ road_gravel_2l : num [1:3100] 0 0 NA 0 0 0 0 0 0 0 ...
## $ truck_trail : num [1:3100] 0.5 0 NA 0 0 0 0 0 0 0 ...
## $ borrowpits : num [1:3100] 0 0 NA 0 0 0 0 0 0 0 ...
## $ sump : num [1:3100] 0 0 NA 0 0 0 0 0 0 0 ...
## $ borrowpit_wet : num [1:3100] 0 0 NA 0 0 0 0 0 0 0 ...
## $ cultivation_abandoned : num [1:3100] 0 0 NA 0 0 0 0 0 0 0 ...
## $ urban_residence : num [1:3100] 0 0 NA 0 0 0 0 0 0 0 ...
## $ country_residence : num [1:3100] 0 0 NA 0 0 0 0 0 0 0 ...
## $ recreation : num [1:3100] 0 0 NA 0 0 0 0 0 0 0 ...
## $ well_other : num [1:3100] 0 0 NA 0 0 0 0 0 0 0 ...
## $ well_bitumen : num [1:3100] 0 0 NA 0 0 0 0 0 0 0 ...
## $ well_cased : num [1:3100] 0 0 NA 0 0 0 0 0 0 0 ...
## $ road_paved_undiv_2l : num [1:3100] 0 0 NA 0 0 0 0 0 0 0 ...
## $ road_unclassified : num [1:3100] 0 0 NA 0 0 0 0 0 0 0 ...
## $ runway : num [1:3100] 0 0 NA 0 0 0 0 0 0 0 ...
## $ clearing_wellpad_unconfirmed: num [1:3100] 0 0 NA 0 0 0 0 0 0 0 ...
## $ facility_unknown : num [1:3100] 0 0 NA 0 0 0 0 0 0 0 ...
## $ borrowpit_dry : num [1:3100] 0 0 NA 0 0 0 0 0 0 0 ...
## $ grvl_sand_pit : num [1:3100] 0 0 NA 0 0 0 0 0 0 0 ...
## $ dugout : num [1:3100] 0 0 NA 0 0 0 0 0 0 0 ...
## $ lagoon : num [1:3100] 0 0 NA 0 0 0 0 0 0 0 ...
## $ open_pit_mine : num [1:3100] 0 0 NA 0 0 0 0 0 0 0 ...
## $ low_impact_seismic : num [1:3100] 0 0 NA 0 0 0 0 0 0 0 ...
## $ surrounding_veg : num [1:3100] 0 0 NA 0 0 0 0 0 0 0 ...
## $ transfer_station : num [1:3100] 0 0 NA 0 0 0 0 0 0 0 ...
## $ facility_other : num [1:3100] 0 0 NA 0 0 0 0 0 0 0 ...
## $ vegetated_edge_railways : num [1:3100] 0 0 NA 0 0 0 0 0 0 0 ...
## $ fruit_vegetables : num [1:3100] 0 0 NA 0 0 0 0 0 0 0 ...
## $ residence_clearing : num [1:3100] 0 0 NA 0 0 0 0 0 0 0 ...
## $ cfo : num [1:3100] 0 0 NA 0 0 0 0 0 0 0 ...
## $ landfill : num [1:3100] 0 0 NA 0 0 0 0 0 0 0 ...
## $ well_cleared_not_confirmed : num [1:3100] 0 0 NA 0 0 0 0 0 0 0 ...
## $ oil_gas_plant : num [1:3100] 0 0 NA 0 0 0 0 0 0 0 ...
## $ urban_industrial : num [1:3100] 0 0 NA 0 0 0 0 0 0 0 ...
## $ road_paved_1l : num [1:3100] 0 0 NA 0 0 0 0 0 0 0 ...
## $ road_paved_undiv_1l : num [1:3100] 0 0 NA 0 0 0 0 0 0 0 ...
## $ road_winter : num [1:3100] 0 0 NA 0 0 0 0 0 0 0 ...
## $ well_cleared_not_drilled : num [1:3100] 0 0 NA 0 0 0 0 0 0 0 ...
## $ well_unknown : num [1:3100] 0 0 NA 0 0 0 0 0 0 0 ...
## $ airp_runway : num [1:3100] 0 0 NA 0 0 0 0 0 0 0 ...
## $ reservoir : num [1:3100] 0 0 NA 0 0 0 0 0 0 0 ...
## $ campground : num [1:3100] 0 0 NA 0 0 0 0 0 0 0 ...
## $ canal : num [1:3100] 0 0 NA 0 0 0 0 0 0 0 ...
## $ camp_industrial : num [1:3100] 0 0 NA 0 0 0 0 0 0 0 ...
## $ rlwy_sgl_track : num [1:3100] 0 0 NA 0 0 0 0 0 0 0 ...
## $ lc_class20 : num [1:3100] 0 0 0.4 0.2 0.5 0 0 0 0 0 ...
## $ lc_class33 : num [1:3100] 0 0 0 0 0 0 0 0 0 0 ...
## $ lc_class34 : num [1:3100] 0 0 0 0 0 0 0 0 0 0 ...
## $ lc_class50 : num [1:3100] 0 0 0.2 0.2 0 0 0 0 0 0 ...
## $ lc_class110 : num [1:3100] 0 0.167 0 0 0 ...
## $ lc_class120 : num [1:3100] 0 0 0 0 0 0 0 0 0 0 ...
## $ lc_class210 : num [1:3100] 0.5 0.333 0.2 0.4 0.5 ...
## $ lc_class220 : num [1:3100] 0.5 0.333 0.2 0.2 0 ...
## $ lc_class230 : num [1:3100] 0 0.167 0 0 0 ...
## - attr(*, "spec")=
## .. cols(
## .. .default = col_number(),
## .. array = col_factor(levels = NULL, ordered = FALSE, include_na = FALSE),
## .. camera = col_factor(levels = NULL, ordered = FALSE, include_na = FALSE),
## .. site = col_factor(levels = NULL, ordered = FALSE, include_na = FALSE),
## .. buff_dist = col_factor(levels = NULL, ordered = FALSE, include_na = FALSE),
## .. harvest_area = col_number(),
## .. crop = col_number(),
## .. well_aband = col_number(),
## .. well_oil = col_number(),
## .. trail = col_number(),
## .. harvest_area_white_zone = col_number(),
## .. conventional_seismic = col_number(),
## .. pipeline = col_number(),
## .. tame_pasture = col_number(),
## .. rough_pasture = col_number(),
## .. rural_residence = col_number(),
## .. transmission_line = col_number(),
## .. well_gas = col_number(),
## .. misc_oil_gas_facility = col_number(),
## .. clearing_unknown = col_number(),
## .. vegetated_edge_roads = col_number(),
## .. road_unimproved = col_number(),
## .. road_gravel_1l = col_number(),
## .. road_gravel_2l = col_number(),
## .. truck_trail = col_number(),
## .. borrowpits = col_number(),
## .. sump = col_number(),
## .. borrowpit_wet = col_number(),
## .. cultivation_abandoned = col_number(),
## .. urban_residence = col_number(),
## .. country_residence = col_number(),
## .. recreation = col_number(),
## .. well_other = col_number(),
## .. well_bitumen = col_number(),
## .. well_cased = col_number(),
## .. road_paved_undiv_2l = col_number(),
## .. road_unclassified = col_number(),
## .. runway = col_number(),
## .. clearing_wellpad_unconfirmed = col_number(),
## .. facility_unknown = col_number(),
## .. borrowpit_dry = col_number(),
## .. grvl_sand_pit = col_number(),
## .. dugout = col_number(),
## .. lagoon = col_number(),
## .. open_pit_mine = col_number(),
## .. low_impact_seismic = col_number(),
## .. surrounding_veg = col_number(),
## .. transfer_station = col_number(),
## .. facility_other = col_number(),
## .. vegetated_edge_railways = col_number(),
## .. fruit_vegetables = col_number(),
## .. residence_clearing = col_number(),
## .. cfo = col_number(),
## .. landfill = col_number(),
## .. well_cleared_not_confirmed = col_number(),
## .. oil_gas_plant = col_number(),
## .. urban_industrial = col_number(),
## .. road_paved_1l = col_number(),
## .. road_paved_undiv_1l = col_number(),
## .. road_winter = col_number(),
## .. well_cleared_not_drilled = col_number(),
## .. well_unknown = col_number(),
## .. airp_runway = col_number(),
## .. reservoir = col_number(),
## .. campground = col_number(),
## .. canal = col_number(),
## .. camp_industrial = col_number(),
## .. rlwy_sgl_track = col_number(),
## .. lc_class20 = col_number(),
## .. lc_class33 = col_number(),
## .. lc_class34 = col_number(),
## .. lc_class50 = col_number(),
## .. lc_class110 = col_number(),
## .. lc_class120 = col_number(),
## .. lc_class210 = col_number(),
## .. lc_class220 = col_number(),
## .. lc_class230 = col_number()
## .. )
## - attr(*, "problems")=<externalptr>Format covariates
There are too many covariates to include in the models individually and many of them describe similar HFI features. We can use the info from the README file in this repository which includes detailed descriptions from the ABMI human footprints wall to wall data download website for Year 2021 OR in the relevant_literature folder of this repository (HFI_2021_v1_0_Metadata_Final.pdf).
the current version of this code for the purposes of the 2022-2023 report used a merged dataset from 2021-2022 and 2022-2023 data, howver each year of data the variables were extracted slightly differenty from GIS so final version of this code will include a different formatting process which will likely occur in the ACME_camera_script_9-2-2024.R or .Rmd
covariates_grouped <- covariates %>%
mutate(borrowpits = rowSums(across(contains('borrowpit'))),
industrial_sites = camp_industrial + oil_gas_plant + open_pit_mine +
rowSums(across(contains('facility'))),
seismic_lines = rowSums(across(contains('seismic'))),
wellsites = rowSums(across(contains('well'))),
roads = rowSums(across(contains('road'))),
havest_areas = rowSums(across(contains('harvest'))),
trails = rowSums(across(contains('trail'))),
residences = rowSums(across(contains('residence'))),
pasture = rowSums(across(contains('pasture'))),
other_transportation_features = runway + airp_runway + rlwy_sgl_track + vegetated_edge_railways,
crops = crop + fruit_vegetables + cultivation_abandoned,
water = lagoon + reservoir + dugout + canal,
.keep = 'unused') %>%
# remove features we don't need
select(!c(recreation,
clearing_unknown,
cfo,
grvl_sand_pit,
transfer_station,
campground,
surrounding_veg,
urban_industrial,
landfill,
sump,
water,
crops,
other_transportation_features,
pasture,
residences
)) %>%
# reorder variables
relocate(c(pipeline,
transmission_line,
borrowpits),
.after = lc_class230)
# see what's left
names(covariates_grouped)## [1] "array" "camera" "site"
## [4] "buff_dist" "lc_class20" "lc_class33"
## [7] "lc_class34" "lc_class50" "lc_class110"
## [10] "lc_class120" "lc_class210" "lc_class220"
## [13] "lc_class230" "pipeline" "transmission_line"
## [16] "borrowpits" "industrial_sites" "seismic_lines"
## [19] "wellsites" "roads" "havest_areas"
## [22] "trails"# check the structure of new data
str(covariates_grouped)## tibble [3,100 × 22] (S3: tbl_df/tbl/data.frame)
## $ array : Factor w/ 4 levels "LU13","LU15",..: 1 1 1 1 1 1 1 1 1 1 ...
## $ camera : Factor w/ 96 levels "18","15","03",..: 1 2 3 4 5 6 7 8 9 10 ...
## $ site : Factor w/ 155 levels "LU13_18","LU13_15",..: 1 2 3 4 5 6 7 8 9 10 ...
## $ buff_dist : Factor w/ 20 levels "250","500","750",..: 1 1 1 1 1 1 1 1 1 1 ...
## $ lc_class20 : num [1:3100] 0 0 0.4 0.2 0.5 0 0 0 0 0 ...
## $ lc_class33 : num [1:3100] 0 0 0 0 0 0 0 0 0 0 ...
## $ lc_class34 : num [1:3100] 0 0 0 0 0 0 0 0 0 0 ...
## $ lc_class50 : num [1:3100] 0 0 0.2 0.2 0 0 0 0 0 0 ...
## $ lc_class110 : num [1:3100] 0 0.167 0 0 0 ...
## $ lc_class120 : num [1:3100] 0 0 0 0 0 0 0 0 0 0 ...
## $ lc_class210 : num [1:3100] 0.5 0.333 0.2 0.4 0.5 ...
## $ lc_class220 : num [1:3100] 0.5 0.333 0.2 0.2 0 ...
## $ lc_class230 : num [1:3100] 0 0.167 0 0 0 ...
## $ pipeline : num [1:3100] 0 0.5 NA 0 0 0 0.5 0 0 0 ...
## $ transmission_line: num [1:3100] 0 0 NA 0 0 0 0 0 0 0 ...
## $ borrowpits : num [1:3100] 0 0 NA 0 0 0 0 0 0 0 ...
## $ industrial_sites : num [1:3100] 0 0 NA 0 0 0 0 0 0 0 ...
## $ seismic_lines : num [1:3100] 0.5 0.5 NA 0.5 1 ...
## $ wellsites : num [1:3100] 0 0 NA 0 0 0 0 0 0 0 ...
## $ roads : num [1:3100] 0 0 NA 0 0 0 0 0 0 0 ...
## $ havest_areas : num [1:3100] 0 0 NA 0 0 0 0 0 0 0 ...
## $ trails : num [1:3100] 0.5 0 NA 0.5 0 ...# check summary of new data
summary(covariates_grouped)## array camera site buff_dist lc_class20
## LU13:820 27 : 80 LU13_18: 20 250 : 155 Min. :0.00000
## LU15:780 32 : 80 LU13_15: 20 500 : 155 1st Qu.:0.00000
## LU21:720 41 : 80 LU13_03: 20 750 : 155 Median :0.05460
## LU01:780 36 : 80 LU13_34: 20 1000 : 155 Mean :0.07304
## 16 : 60 LU13_57: 20 1250 : 155 3rd Qu.:0.11321
## 21 : 60 LU13_16: 20 1500 : 155 Max. :0.59091
## (Other):2660 (Other):2980 (Other):2170
## lc_class33 lc_class34 lc_class50 lc_class110
## Min. :0.0000000 Min. :0.00000 Min. :0.0000 Min. :0.00000
## 1st Qu.:0.0000000 1st Qu.:0.00000 1st Qu.:0.1000 1st Qu.:0.05882
## Median :0.0000000 Median :0.00000 Median :0.1686 Median :0.12500
## Mean :0.0005196 Mean :0.01445 Mean :0.1829 Mean :0.12674
## 3rd Qu.:0.0000000 3rd Qu.:0.01316 3rd Qu.:0.2500 3rd Qu.:0.17550
## Max. :0.0909091 Max. :0.33333 Max. :1.0000 Max. :0.55556
##
## lc_class120 lc_class210 lc_class220 lc_class230
## Min. :0.00000 Min. :0.0000 Min. :0.0000 Min. :0.00000
## 1st Qu.:0.00000 1st Qu.:0.1064 1st Qu.:0.1726 1st Qu.:0.08217
## Median :0.00000 Median :0.1591 Median :0.2222 Median :0.15152
## Mean :0.03835 Mean :0.1813 Mean :0.2268 Mean :0.15583
## 3rd Qu.:0.00000 3rd Qu.:0.2179 3rd Qu.:0.2735 3rd Qu.:0.22093
## Max. :1.00000 Max. :1.0000 Max. :1.0000 Max. :0.66667
##
## pipeline transmission_line borrowpits industrial_sites
## Min. :0.00000 Min. :0.000000 Min. :0.00000 Min. :0.000000
## 1st Qu.:0.00000 1st Qu.:0.000000 1st Qu.:0.00000 1st Qu.:0.000000
## Median :0.03637 Median :0.000000 Median :0.00000 Median :0.000000
## Mean :0.06940 Mean :0.004712 Mean :0.01108 Mean :0.001448
## 3rd Qu.:0.10638 3rd Qu.:0.000000 3rd Qu.:0.01807 3rd Qu.:0.000000
## Max. :1.00000 Max. :0.500000 Max. :0.16667 Max. :0.111111
## NA's :8 NA's :8 NA's :8 NA's :8
## seismic_lines wellsites roads havest_areas
## Min. :0.0000 Min. :0.00000 Min. :0.00000 Min. :0.00000
## 1st Qu.:0.2774 1st Qu.:0.01541 1st Qu.:0.00000 1st Qu.:0.00000
## Median :0.3868 Median :0.04408 Median :0.05939 Median :0.00000
## Mean :0.4173 Mean :0.05748 Mean :0.15189 Mean :0.04801
## 3rd Qu.:0.5000 3rd Qu.:0.08125 3rd Qu.:0.27978 3rd Qu.:0.04506
## Max. :1.0000 Max. :0.50000 Max. :0.83333 Max. :0.83333
## NA's :8 NA's :8 NA's :8 NA's :8
## trails
## Min. :0.0000
## 1st Qu.:0.0617
## Median :0.1522
## Mean :0.1874
## 3rd Qu.:0.2712
## Max. :1.0000
## NA's :8# there are some NAs in the data which will cause problems with modeling/visualization of data ignore for now but will explore these sites specifically after report
covariates_grouped <- covariates_grouped %>%
# remove rows with NAs
na.omit()Subset data & correlation plots
Marissa try to get the purrr code for this to work later
Now we need to subset the data for each buffer width, and then in the same loop let’s make correlation plots for these variables within each buffer
# Couldn't get this to work in purrr yet so using a loop to subset the data, create the plots, and save them all in one section... NEAT
buffer_frames<-list()
for (i in unique(covariates_grouped$buff_dist)){
print(i)
#Subset data based on radius
df<-covariates_grouped%>%
filter(buff_dist == i)
#rename dataframe on the fly
assign(paste("df", i, sep ="_"), df)
#list of dataframes
buffer_frames<-c(buffer_frames, list(df))
#Subset data based on radius
df<-covariates_grouped%>%
filter(buff_dist == i)%>%
select(where(is.numeric))
#compute a correlation matrix (watch for errors)
matrix<-cor(df)
#print and save the correlation plot on the go
#renaming for each buffer as we do
png(file.path("figures/", paste("correlation_", i, ".png")))
corrplot::corrplot(matrix,
type = 'upper',
tl.col = 'black',
title = paste0('Variable correlation plot at ', i))
dev.off()
}## [1] "250"## Warning in cor(df): the standard deviation is zero## [1] "500"## Warning in cor(df): the standard deviation is zero## [1] "750"## [1] "1000"## [1] "1250"## [1] "1500"## [1] "1750"## [1] "2000"## [1] "2250"## [1] "2500"## [1] "2750"## [1] "3000"## [1] "3250"## [1] "3500"## [1] "3750"## [1] "4000"## [1] "4250"## [1] "4500"## [1] "4750"## [1] "5000"# name list objects so we can extract names for plotting
buffer_frames <- buffer_frames %>%
# absurdly long way to do this but for sake of time fuck it
purrr::set_names('250 meter buffer',
'500 meter buffer',
'750 meter buffer',
'1000 meter buffer',
'1250 meter buffer',
'1500 meter buffer',
'1750 meter buffer',
'2000 meter buffer',
'2250 meter buffer',
'2500 meter buffer',
'2750 meter buffer',
'3000 meter buffer',
'3250 meter buffer',
'3500 meter buffer',
'3750 meter buffer',
'4000 meter buffer',
'4250 meter buffer',
'4500 meter buffer',
'4750 meter buffer',
'5000 meter buffer')Exploratory plots
add more to this section in later when we have more time to explore the covariates and choose which should be inlcuded etc.
hfi_histograms <- buffer_frames %>%
purrr::imap(
~.x %>%
# filter to just the HFI variables
select(where(is.numeric) &
! starts_with('lc_class')) %>%
# pipe into hist.data.frame function to make histograms for each variable
hist.data.frame(mtitl = paste0('Histograms of HFI variables at ', .y)))
Now let’s do the same thing with the landcover variables
lc_histograms <- buffer_frames %>%
purrr::imap(
~.x %>%
# filter to just the landcover variables
select(where(is.numeric) &
starts_with('lc_class')) %>%
# pipe into hist.data.frame function to make histograms for each variable
hist.data.frame(mtitl = paste0('Histograms of landcover variables at ', .y)))
Add response metric
Now that we have the covariate data formatted we need to add the response metric (monthly proportional presence/absence) to the data frames
final_df <- buffer_frames %>%
purrr::map(
~.x %>%
left_join(prop_detections,
by = 'site'))Analysis
Now we are going to run a global model which includes all HFI and LC variables that at first glance (will do a more thorough check later) seem to have enough data to include as covariates for each buffer width, and then we will compare these models see which buffer width best fit the data for each species.
We don’t need to do ALL the species since many don’t have enough data.
Refer to the ACME_camera_script_9-2-2024.html or .Rmd the plot for proportional monthly detections should provide info on which species we have enough data for, can be found under Response metrics/3.Proportion monthly detections
A brief look at this fig indicates that we have enough for all the mammals in the prop_detections data frame except
- cougar
- wolverine
- caribou??? (may have enough may not)
Black bear
ANDREW/MARISSA FIX later
there is probably a way to shorten the following code to select particular species, I saw Andrew’s for loop in the draft script he wrote but couldn’t quite figure it out so I did this instead, maybe we can merge approaches?
Global model
Let’s start with bears and use purrr to create a global model for every buffer distance
black_bear_mods <- final_df %>%
# use purrr map to fun the same functions for all buffer sizes ((all objects in list))
purrr::map(
~.x %>%
# glmmTMB function let's us run the proportional binomial model using cbind to combine the present and absent columns for each species
glmmTMB::glmmTMB(cbind(black_bear, absent_black_bear) ~
seismic_lines +
pipeline +
borrowpits +
wellsites +
roads +
trails +
lc_class20 +
lc_class34 +
lc_class50 +
lc_class110 +
lc_class210 +
lc_class220 +
lc_class230 +
(1|array),
data = .,
family = 'binomial'))Model selection
We will use the model.sel() function from the
MuMIn package to compare the global models for each buffer
width and see which buffer fits the bear data best
model.sel(black_bear_mods)## Warning in model.sel.default(black_bear_mods): models are not all fitted to the
## same data## Model selection table
## cnd((Int)) dsp((Int)) cnd(brr) cnd(lc_c11) cnd(lc_c20)
## 250 meter buffer -0.8237 + 1.0280 0.40780 0.5002000
## 500 meter buffer -1.8100 + -8.8800 1.48300 2.1650000
## 1000 meter buffer -1.2440 + 5.0040 -0.79470 0.9692000
## 750 meter buffer -0.9882 + 3.2100 0.08085 1.0990000
## 1500 meter buffer -1.2430 + -1.6840 -1.97700 -0.1470000
## 4000 meter buffer -2.5810 + -10.0700 1.06300 2.4720000
## 4250 meter buffer -3.1640 + -12.3600 0.46190 3.0210000
## 1250 meter buffer -1.5400 + -3.2830 -0.87330 0.6002000
## 4500 meter buffer -3.1770 + -11.5900 -0.12570 2.3210000
## 5000 meter buffer -2.9500 + -9.5120 2.01100 3.0470000
## 4750 meter buffer -2.9020 + -10.2100 1.21200 2.9250000
## 3750 meter buffer -2.0630 + -7.5570 -0.28530 1.8680000
## 1750 meter buffer -1.2670 + -2.4830 -1.42300 0.1190000
## 3500 meter buffer -2.1010 + -7.2870 -0.70810 1.3820000
## 3250 meter buffer -0.3758 + -3.9270 -2.37000 0.0878400
## 3000 meter buffer -0.2031 + -3.0550 -2.93400 -0.1205000
## 2750 meter buffer 0.2775 + 0.8066 -2.37500 -0.2794000
## 2250 meter buffer -0.6461 + 2.6130 -1.06400 0.5436000
## 2000 meter buffer -0.7442 + 0.5012 -1.60200 -0.0001237
## 2500 meter buffer -0.1102 + 0.9909 -2.07700 -0.2207000
## cnd(lc_c21) cnd(lc_c22) cnd(lc_c23) cnd(lc_c34) cnd(lc_c50)
## 250 meter buffer -0.02139 0.82310 0.334200 -0.01718 0.05003
## 500 meter buffer 1.02300 2.56100 1.662000 2.70000 1.47400
## 1000 meter buffer -0.99610 0.63190 0.488000 -0.19810 0.25530
## 750 meter buffer -0.68520 1.29900 0.121600 -1.50100 0.01754
## 1500 meter buffer -1.13400 1.20900 -0.657200 -2.59900 -0.55680
## 4000 meter buffer -1.28200 1.08200 3.415000 31.62000 1.35600
## 4250 meter buffer -1.30100 2.05700 3.583000 26.01000 1.64500
## 1250 meter buffer -0.56750 1.28500 0.170100 0.45260 0.06635
## 4500 meter buffer -2.24200 1.93700 3.601000 -1.36700 1.68800
## 5000 meter buffer -1.26500 0.35690 5.558000 24.91000 2.14500
## 4750 meter buffer -1.73800 0.94160 4.978000 22.74000 1.84500
## 3750 meter buffer -1.88700 0.39800 2.642000 11.30000 1.01100
## 1750 meter buffer -1.02700 1.33500 0.063110 4.37900 0.10100
## 3500 meter buffer -2.00500 0.60230 2.294000 5.30400 0.86070
## 3250 meter buffer -2.44600 -0.85350 -0.697700 -4.48900 -0.64980
## 3000 meter buffer -2.13700 -0.80130 -1.454000 -8.13900 -1.22600
## 2750 meter buffer -1.95800 -1.51600 -1.743000 -10.23000 -1.66000
## 2250 meter buffer -1.75100 -0.54820 0.076010 0.62950 -0.64770
## 2000 meter buffer -1.28200 0.01076 -0.002506 0.43760 -0.64780
## 2500 meter buffer -1.49300 -1.10300 -1.262000 0.99500 -1.82600
## cnd(ppl) cnd(rds) cnd(ssm_lns) cnd(trl) cnd(wll) df logLik
## 250 meter buffer -0.58350 0.04315 -0.24170 -0.1009 -1.53600 15 -294.398
## 500 meter buffer -1.18100 -0.37350 -0.52020 -0.5640 0.49210 15 -313.151
## 1000 meter buffer 0.46020 0.66440 0.71130 0.6546 -1.01000 15 -313.579
## 750 meter buffer -0.33060 0.39990 0.08933 0.1862 -0.04885 15 -313.611
## 1500 meter buffer 1.03500 0.96360 1.27100 1.6530 -0.36790 15 -314.163
## 4000 meter buffer -0.05392 0.81850 0.80620 1.1670 -0.69150 15 -315.277
## 4250 meter buffer 0.33400 1.50500 1.15300 1.4560 0.31490 15 -315.758
## 1250 meter buffer 0.95440 0.90300 0.84580 1.2300 -0.93390 15 -315.859
## 4500 meter buffer 1.63800 2.26200 1.74200 1.8280 -0.70290 15 -316.302
## 5000 meter buffer -1.26300 1.42800 0.56130 0.5805 0.78170 15 -317.172
## 4750 meter buffer -0.63540 1.47000 0.79400 0.7656 0.83780 15 -317.404
## 3750 meter buffer 0.57060 1.07800 1.07100 1.3880 -0.48320 15 -317.432
## 1750 meter buffer 0.40550 0.03219 0.92360 0.6467 1.10800 15 -317.524
## 3500 meter buffer 0.72160 1.32200 1.38500 1.4750 0.14660 15 -318.383
## 3250 meter buffer 0.18050 0.74830 1.18900 1.1820 1.57800 15 -318.817
## 3000 meter buffer 0.46020 0.92300 1.24900 1.3140 2.00700 15 -319.422
## 2750 meter buffer 0.87090 0.47250 0.88000 0.8917 1.42000 15 -320.792
## 2250 meter buffer 0.99960 0.23030 0.74480 0.8051 0.13300 15 -321.218
## 2000 meter buffer 0.76620 0.27650 0.88490 0.7060 0.53110 15 -321.295
## 2500 meter buffer 0.59280 0.25370 0.97210 1.0170 1.66700 15 -321.567
## AICc delta weight
## 250 meter buffer 622.5 0.00 1
## 500 meter buffer 659.8 37.29 0
## 1000 meter buffer 660.7 38.14 0
## 750 meter buffer 660.8 38.21 0
## 1500 meter buffer 661.9 39.31 0
## 4000 meter buffer 664.1 41.54 0
## 4250 meter buffer 665.0 42.50 0
## 1250 meter buffer 665.2 42.70 0
## 4500 meter buffer 666.1 43.59 0
## 5000 meter buffer 667.9 45.33 0
## 4750 meter buffer 668.3 45.79 0
## 3750 meter buffer 668.4 45.85 0
## 1750 meter buffer 668.6 46.03 0
## 3500 meter buffer 670.3 47.75 0
## 3250 meter buffer 671.2 48.62 0
## 3000 meter buffer 672.4 49.83 0
## 2750 meter buffer 675.1 52.57 0
## 2250 meter buffer 676.0 53.42 0
## 2000 meter buffer 676.1 53.58 0
## 2500 meter buffer 676.7 54.12 0
## Models ranked by AICc(x)
## Random terms (all models):
## cond(1 | array)hmmmm seems fishy to me that the 250 meter buffer which is the only one that had missing data would perform THAT much better than all the others, and really you shouldn’t compare models if they aren’t run on the same data, hence the warning message
Let’s remove the 250 buffer and see what happens
black_bear_mods_no250 <- black_bear_mods %>%
# purrr::discard_at will remove an item from a list
purrr::discard_at('250 meter buffer')
# run model selection again
model.sel(black_bear_mods_no250)## Model selection table
## cnd((Int)) dsp((Int)) cnd(brr) cnd(lc_c11) cnd(lc_c20)
## 500 meter buffer -1.8100 + -8.8800 1.48300 2.1650000
## 1000 meter buffer -1.2440 + 5.0040 -0.79470 0.9692000
## 750 meter buffer -0.9882 + 3.2100 0.08085 1.0990000
## 1500 meter buffer -1.2430 + -1.6840 -1.97700 -0.1470000
## 4000 meter buffer -2.5810 + -10.0700 1.06300 2.4720000
## 4250 meter buffer -3.1640 + -12.3600 0.46190 3.0210000
## 1250 meter buffer -1.5400 + -3.2830 -0.87330 0.6002000
## 4500 meter buffer -3.1770 + -11.5900 -0.12570 2.3210000
## 5000 meter buffer -2.9500 + -9.5120 2.01100 3.0470000
## 4750 meter buffer -2.9020 + -10.2100 1.21200 2.9250000
## 3750 meter buffer -2.0630 + -7.5570 -0.28530 1.8680000
## 1750 meter buffer -1.2670 + -2.4830 -1.42300 0.1190000
## 3500 meter buffer -2.1010 + -7.2870 -0.70810 1.3820000
## 3250 meter buffer -0.3758 + -3.9270 -2.37000 0.0878400
## 3000 meter buffer -0.2031 + -3.0550 -2.93400 -0.1205000
## 2750 meter buffer 0.2775 + 0.8066 -2.37500 -0.2794000
## 2250 meter buffer -0.6461 + 2.6130 -1.06400 0.5436000
## 2000 meter buffer -0.7442 + 0.5012 -1.60200 -0.0001237
## 2500 meter buffer -0.1102 + 0.9909 -2.07700 -0.2207000
## cnd(lc_c21) cnd(lc_c22) cnd(lc_c23) cnd(lc_c34) cnd(lc_c50)
## 500 meter buffer 1.0230 2.56100 1.662000 2.7000 1.47400
## 1000 meter buffer -0.9961 0.63190 0.488000 -0.1981 0.25530
## 750 meter buffer -0.6852 1.29900 0.121600 -1.5010 0.01754
## 1500 meter buffer -1.1340 1.20900 -0.657200 -2.5990 -0.55680
## 4000 meter buffer -1.2820 1.08200 3.415000 31.6200 1.35600
## 4250 meter buffer -1.3010 2.05700 3.583000 26.0100 1.64500
## 1250 meter buffer -0.5675 1.28500 0.170100 0.4526 0.06635
## 4500 meter buffer -2.2420 1.93700 3.601000 -1.3670 1.68800
## 5000 meter buffer -1.2650 0.35690 5.558000 24.9100 2.14500
## 4750 meter buffer -1.7380 0.94160 4.978000 22.7400 1.84500
## 3750 meter buffer -1.8870 0.39800 2.642000 11.3000 1.01100
## 1750 meter buffer -1.0270 1.33500 0.063110 4.3790 0.10100
## 3500 meter buffer -2.0050 0.60230 2.294000 5.3040 0.86070
## 3250 meter buffer -2.4460 -0.85350 -0.697700 -4.4890 -0.64980
## 3000 meter buffer -2.1370 -0.80130 -1.454000 -8.1390 -1.22600
## 2750 meter buffer -1.9580 -1.51600 -1.743000 -10.2300 -1.66000
## 2250 meter buffer -1.7510 -0.54820 0.076010 0.6295 -0.64770
## 2000 meter buffer -1.2820 0.01076 -0.002506 0.4376 -0.64780
## 2500 meter buffer -1.4930 -1.10300 -1.262000 0.9950 -1.82600
## cnd(ppl) cnd(rds) cnd(ssm_lns) cnd(trl) cnd(wll) df logLik
## 500 meter buffer -1.18100 -0.37350 -0.52020 -0.5640 0.49210 15 -313.151
## 1000 meter buffer 0.46020 0.66440 0.71130 0.6546 -1.01000 15 -313.579
## 750 meter buffer -0.33060 0.39990 0.08933 0.1862 -0.04885 15 -313.611
## 1500 meter buffer 1.03500 0.96360 1.27100 1.6530 -0.36790 15 -314.163
## 4000 meter buffer -0.05392 0.81850 0.80620 1.1670 -0.69150 15 -315.277
## 4250 meter buffer 0.33400 1.50500 1.15300 1.4560 0.31490 15 -315.758
## 1250 meter buffer 0.95440 0.90300 0.84580 1.2300 -0.93390 15 -315.859
## 4500 meter buffer 1.63800 2.26200 1.74200 1.8280 -0.70290 15 -316.302
## 5000 meter buffer -1.26300 1.42800 0.56130 0.5805 0.78170 15 -317.172
## 4750 meter buffer -0.63540 1.47000 0.79400 0.7656 0.83780 15 -317.404
## 3750 meter buffer 0.57060 1.07800 1.07100 1.3880 -0.48320 15 -317.432
## 1750 meter buffer 0.40550 0.03219 0.92360 0.6467 1.10800 15 -317.524
## 3500 meter buffer 0.72160 1.32200 1.38500 1.4750 0.14660 15 -318.383
## 3250 meter buffer 0.18050 0.74830 1.18900 1.1820 1.57800 15 -318.817
## 3000 meter buffer 0.46020 0.92300 1.24900 1.3140 2.00700 15 -319.422
## 2750 meter buffer 0.87090 0.47250 0.88000 0.8917 1.42000 15 -320.792
## 2250 meter buffer 0.99960 0.23030 0.74480 0.8051 0.13300 15 -321.218
## 2000 meter buffer 0.76620 0.27650 0.88490 0.7060 0.53110 15 -321.295
## 2500 meter buffer 0.59280 0.25370 0.97210 1.0170 1.66700 15 -321.567
## AICc delta weight
## 500 meter buffer 659.8 0.00 0.331
## 1000 meter buffer 660.7 0.86 0.216
## 750 meter buffer 660.8 0.92 0.209
## 1500 meter buffer 661.9 2.02 0.120
## 4000 meter buffer 664.1 4.25 0.040
## 4250 meter buffer 665.0 5.21 0.024
## 1250 meter buffer 665.2 5.42 0.022
## 4500 meter buffer 666.1 6.30 0.014
## 5000 meter buffer 667.9 8.04 0.006
## 4750 meter buffer 668.3 8.51 0.005
## 3750 meter buffer 668.4 8.56 0.005
## 1750 meter buffer 668.6 8.75 0.004
## 3500 meter buffer 670.3 10.46 0.002
## 3250 meter buffer 671.2 11.33 0.001
## 3000 meter buffer 672.4 12.54 0.001
## 2750 meter buffer 675.1 15.28 0.000
## 2250 meter buffer 676.0 16.13 0.000
## 2000 meter buffer 676.1 16.29 0.000
## 2500 meter buffer 676.7 16.83 0.000
## Models ranked by AICc(x)
## Random terms (all models):
## cond(1 | array)this looks much more realistic; the 500 m buffer is top model for black bears
So what we will do for each species is remove the 250 meter buffer for now since there are some data missing. and compare just the other buffer sizes that contain the full data set
Model summary
Let’s take a look at the model summary for the top model
summary(black_bear_mods_no250$`500 meter buffer`)## Family: binomial ( logit )
## Formula:
## cbind(black_bear, absent_black_bear) ~ seismic_lines + pipeline +
## borrowpits + wellsites + roads + trails + lc_class20 + lc_class34 +
## lc_class50 + lc_class110 + lc_class210 + lc_class220 + lc_class230 +
## (1 | array)
## Data: .
##
## AIC BIC logLik deviance df.resid
## 656.3 701.7 -313.2 626.3 137
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## array (Intercept) 0.02745 0.1657
## Number of obs: 152, groups: array, 4
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.8105 0.6513 -2.780 0.00544 **
## seismic_lines -0.5202 0.4415 -1.178 0.23865
## pipeline -1.1814 0.6352 -1.860 0.06292 .
## borrowpits -8.8805 3.8839 -2.287 0.02222 *
## wellsites 0.4921 0.8628 0.570 0.56845
## roads -0.3735 0.5927 -0.630 0.52862
## trails -0.5640 0.4732 -1.192 0.23325
## lc_class20 2.1652 0.9060 2.390 0.01686 *
## lc_class34 2.7001 1.4844 1.819 0.06890 .
## lc_class50 1.4742 0.7802 1.890 0.05883 .
## lc_class110 1.4826 0.8658 1.712 0.08684 .
## lc_class210 1.0230 0.7283 1.405 0.16014
## lc_class220 2.5607 0.8105 3.160 0.00158 **
## lc_class230 1.6625 0.8342 1.993 0.04627 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Let’s repeat this process for each species that we have enough data for.
Caribou
We may or may not have enough data for caribou but let’s try it at least for this preliminary report
We can use the same code from black beasr (above) to run global models for each buffer width except remember we want to remove 250 meters
And in the same chunk to save time let’s also run the
model.sel() function
caribou_mods_no250 <- final_df %>%
# remove 350 meter buffer width
purrr::discard_at('250 meter buffer') %>%
# use purrr map to make global models for all other buffer sizes
purrr::map(
~.x %>%
glmmTMB::glmmTMB(cbind(caribou, absent_caribou) ~
seismic_lines +
pipeline +
borrowpits +
wellsites +
roads +
trails +
lc_class20 +
lc_class34 +
lc_class50 +
lc_class110 +
lc_class210 +
lc_class220 +
lc_class230 +
(1|array),
data = .,
family = 'binomial'))## Warning in finalizeTMB(TMBStruc, obj, fit, h, data.tmb.old): Model convergence
## problem; singular convergence (7). See vignette('troubleshooting'),
## help('diagnose')
## Warning in finalizeTMB(TMBStruc, obj, fit, h, data.tmb.old): Model convergence
## problem; singular convergence (7). See vignette('troubleshooting'),
## help('diagnose')# model selection
model.sel(caribou_mods_no250)## Model selection table
## cnd((Int)) dsp((Int)) cnd(brr) cnd(lc_c11) cnd(lc_c20)
## 1250 meter buffer -14.540 + 1.5550 1.813 2.7510
## 2250 meter buffer -14.910 + -42.4600 7.942 5.9510
## 750 meter buffer -11.580 + 14.0000 1.297 0.8007
## 4750 meter buffer -45.610 + -78.1400 47.080 20.1100
## 4500 meter buffer -45.040 + -50.1400 30.330 16.0500
## 1000 meter buffer -17.550 + 4.9460 5.340 5.1840
## 2500 meter buffer -18.790 + -58.7900 13.190 7.3110
## 2000 meter buffer -9.641 + -14.6300 4.367 1.0910
## 3750 meter buffer -15.180 + -4.3390 11.910 -0.7118
## 5000 meter buffer -33.410 + -73.3100 40.230 9.5270
## 3250 meter buffer -8.524 + 31.6800 11.680 -3.6230
## 1750 meter buffer -12.710 + 0.5022 6.596 7.0970
## 4250 meter buffer -28.480 + -32.9700 22.640 11.0600
## 1500 meter buffer -14.540 + -1.2960 4.483 4.7300
## 3500 meter buffer -13.880 + 15.0600 15.540 4.6380
## 4000 meter buffer -23.140 + -31.0600 21.470 9.5710
## 500 meter buffer -10.320 + 7.9850 5.291 7.4950
## 2750 meter buffer -12.800 + -7.7580 9.312 4.9720
## 3000 meter buffer -11.720 + 8.4750 10.870 5.3860
## cnd(lc_c21) cnd(lc_c22) cnd(lc_c23) cnd(lc_c34) cnd(lc_c50)
## 1250 meter buffer 4.1770 9.148 10.180 -17.300 7.456
## 2250 meter buffer 11.6000 9.277 6.147 58.240 15.970
## 750 meter buffer 0.9352 4.131 5.492 5.269 1.368
## 4750 meter buffer 15.6100 32.750 45.230 236.000 35.520
## 4500 meter buffer 13.3300 37.440 20.770 160.800 29.410
## 1000 meter buffer 7.2680 12.310 10.310 -3.551 8.005
## 2500 meter buffer 12.4600 13.760 11.000 75.300 19.860
## 2000 meter buffer 4.9070 5.805 2.698 2.155 10.590
## 3750 meter buffer 1.8400 11.600 15.950 13.000 10.500
## 5000 meter buffer 3.9200 21.730 35.550 449.900 24.070
## 3250 meter buffer 0.8793 1.149 6.083 -17.350 3.152
## 1750 meter buffer 9.6800 8.712 6.713 4.289 10.260
## 4250 meter buffer 11.0300 25.120 17.630 111.300 20.290
## 1500 meter buffer 5.4220 7.313 12.420 20.860 8.929
## 3500 meter buffer 4.8220 11.840 14.260 40.310 8.349
## 4000 meter buffer 8.5250 21.090 19.330 95.270 19.250
## 500 meter buffer 3.5200 6.002 6.685 5.520 5.465
## 2750 meter buffer 6.2850 4.502 10.210 26.040 12.090
## 3000 meter buffer 6.2060 3.140 11.420 35.100 9.028
## cnd(ppl) cnd(rds) cnd(ssm_lns) cnd(trl) cnd(wll) df logLik
## 1250 meter buffer -2.218e+01 6.0410 1.21100 7.779 2.068 15 -78.933
## 2250 meter buffer 6.762e+00 -1.6090 -2.36000 6.273 -23.730 15 -80.710
## 750 meter buffer -6.771e+02 4.1230 3.94700 6.967 -0.580 15 -81.359
## 4750 meter buffer -6.312e+00 10.6400 5.41500 27.770 -39.350 15 -81.431
## 4500 meter buffer -1.013e+01 19.2700 13.66000 31.880 -11.500 15 -81.905
## 1000 meter buffer -1.391e+01 6.8920 2.33800 8.140 3.393 15 -82.021
## 2500 meter buffer 1.456e+00 -0.1749 -1.71000 5.894 -13.540 15 -82.836
## 2000 meter buffer 6.652e-01 -3.1290 -2.23300 5.616 -18.860 15 -83.341
## 3750 meter buffer -5.895e+00 -2.8270 0.61120 8.805 -14.770 15 -84.593
## 5000 meter buffer 1.521e+00 1.5530 0.67270 24.390 -36.590 15 -84.700
## 3250 meter buffer -9.612e+00 -3.7050 -0.28470 5.109 -6.538 15 -84.884
## 1750 meter buffer -1.051e+01 0.2878 -2.30000 5.400 -10.960 15 -85.338
## 4250 meter buffer -8.175e+00 5.4900 3.53400 16.660 -12.770 15 -85.748
## 1500 meter buffer -1.330e+01 1.7300 -0.74690 5.934 5.523 15 -85.796
## 3500 meter buffer -1.238e+01 -4.3080 -2.43600 3.322 -3.128 15 -85.798
## 4000 meter buffer -8.715e+00 0.3123 -0.59530 10.270 -12.020 15 -86.346
## 500 meter buffer -1.070e+03 1.7350 0.07777 1.639 -2.600 15 -87.467
## 2750 meter buffer -5.968e-02 -1.1960 -1.63800 5.990 -16.510 15 -87.585
## 3000 meter buffer -7.765e+00 -1.8140 -2.30000 3.983 -8.027 15 -89.041
## AICc delta weight
## 1250 meter buffer 191.4 0.00 0.674
## 2250 meter buffer 194.9 3.55 0.114
## 750 meter buffer 196.2 4.85 0.060
## 4750 meter buffer 196.4 4.99 0.055
## 4500 meter buffer 197.3 5.94 0.035
## 1000 meter buffer 197.6 6.18 0.031
## 2500 meter buffer 199.2 7.81 0.014
## 2000 meter buffer 200.2 8.82 0.008
## 3750 meter buffer 202.7 11.32 0.002
## 5000 meter buffer 202.9 11.53 0.002
## 3250 meter buffer 203.3 11.90 0.002
## 1750 meter buffer 204.2 12.81 0.001
## 4250 meter buffer 205.0 13.63 0.001
## 1500 meter buffer 205.1 13.72 0.001
## 3500 meter buffer 205.1 13.73 0.001
## 4000 meter buffer 206.2 14.82 0.000
## 500 meter buffer 208.5 17.07 0.000
## 2750 meter buffer 208.7 17.30 0.000
## 3000 meter buffer 211.6 20.22 0.000
## Models ranked by AICc(x)
## Random terms (all models):
## cond(1 | array)We get a warning that there are some model convergence problems, I expect this is because we don’t have enough data for caribou but I don’t have time to dig into this now so we will investigate more closely for final analysis
For caribou 1250m buffer is top model for now
Let’s take a closer look at the top model summary
summary(caribou_mods_no250$`1250 meter buffer`)## Family: binomial ( logit )
## Formula:
## cbind(caribou, absent_caribou) ~ seismic_lines + pipeline + borrowpits +
## wellsites + roads + trails + lc_class20 + lc_class34 + lc_class50 +
## lc_class110 + lc_class210 + lc_class220 + lc_class230 + (1 | array)
## Data: .
##
## AIC BIC logLik deviance df.resid
## 187.9 233.2 -78.9 157.9 137
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## array (Intercept) 11.55 3.398
## Number of obs: 152, groups: array, 4
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -14.542 3.972 -3.661 0.000251 ***
## seismic_lines 1.211 2.657 0.456 0.648624
## pipeline -22.176 8.375 -2.648 0.008096 **
## borrowpits 1.555 12.097 0.129 0.897737
## wellsites 2.068 5.003 0.413 0.679401
## roads 6.041 2.905 2.080 0.037535 *
## trails 7.779 3.021 2.575 0.010027 *
## lc_class20 2.751 4.320 0.637 0.524302
## lc_class34 -17.298 14.425 -1.199 0.230453
## lc_class50 7.456 4.250 1.754 0.079403 .
## lc_class110 1.813 4.194 0.432 0.665533
## lc_class210 4.177 4.049 1.032 0.302239
## lc_class220 9.148 4.099 2.232 0.025645 *
## lc_class230 10.178 4.669 2.180 0.029270 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1There’s nothing that catches my eye immediately as being sus about this particular model so it may not have been one with convergence issues. We will keep it in report for now
Coyote
coyote_mods_no250 <- final_df %>%
# remove 350 meter buffer width
purrr::discard_at('250 meter buffer') %>%
# use purrr map to make global models for all other buffer sizes
purrr::map(
~.x %>%
glmmTMB::glmmTMB(cbind(coyote, absent_coyote) ~
seismic_lines +
pipeline +
borrowpits +
wellsites +
roads +
trails +
lc_class20 +
lc_class34 +
lc_class50 +
lc_class110 +
lc_class210 +
lc_class220 +
lc_class230 +
(1|array),
data = .,
family = 'binomial'))
# model selection
model.sel(coyote_mods_no250)## Model selection table
## cnd((Int)) dsp((Int)) cnd(brr) cnd(lc_c11) cnd(lc_c20)
## 4500 meter buffer -0.60210 + -2.9490 -4.1190 -1.2980
## 4750 meter buffer 0.03242 + -6.6400 -6.0680 -0.6162
## 5000 meter buffer 0.83790 + -4.0080 -5.7300 -0.2184
## 4250 meter buffer -0.14980 + -3.3160 -3.9370 -1.2580
## 3750 meter buffer -0.57010 + -4.9410 -6.0980 -1.6480
## 3500 meter buffer -0.60770 + -5.7420 -5.3230 -2.9050
## 4000 meter buffer -1.50200 + -7.1110 -3.8650 -0.7927
## 3000 meter buffer 0.43020 + -7.6520 -6.2080 -3.3750
## 2750 meter buffer 0.79040 + -6.1840 -6.7980 -4.0940
## 3250 meter buffer 0.04700 + -8.4530 -5.9800 -2.8080
## 2500 meter buffer -0.88790 + -5.6500 -4.6990 -2.3540
## 2250 meter buffer -0.76660 + -2.6300 -4.3370 -1.9770
## 2000 meter buffer -0.88490 + -1.2490 -4.0270 -2.3010
## 1750 meter buffer -1.65100 + -4.7080 -2.5830 -1.0160
## 1500 meter buffer -2.20600 + -5.8600 -1.9320 -1.5950
## 500 meter buffer -2.11300 + -8.9480 -0.3059 2.3440
## 1250 meter buffer -1.61300 + -2.9500 -2.7180 -1.5740
## 1000 meter buffer -1.77500 + -0.7574 -1.0300 -0.3119
## 750 meter buffer -1.18000 + 1.6630 -1.2320 0.1889
## cnd(lc_c21) cnd(lc_c22) cnd(lc_c23) cnd(lc_c34) cnd(lc_c50)
## 4500 meter buffer -4.57800 -4.2840 4.2890 -80.3300 -5.674000
## 4750 meter buffer -4.70100 -4.4020 4.0110 -74.9200 -5.951000
## 5000 meter buffer -3.93100 -5.7880 3.7670 -76.5600 -5.802000
## 4250 meter buffer -3.37000 -4.0610 4.1160 -31.3800 -4.728000
## 3750 meter buffer -4.62800 -3.8570 2.9060 -40.4000 -5.039000
## 3500 meter buffer -4.21100 -3.2990 3.0520 -31.9400 -3.942000
## 4000 meter buffer -3.70200 -2.1620 4.5870 -23.6500 -4.788000
## 3000 meter buffer -3.24600 -3.4210 0.2559 -16.4200 -4.045000
## 2750 meter buffer -2.55400 -3.9600 -1.1930 -23.6900 -3.853000
## 3250 meter buffer -3.20400 -3.3130 1.4690 -18.9000 -3.317000
## 2500 meter buffer -1.78700 -2.5510 0.8870 -14.1700 -2.907000
## 2250 meter buffer -1.54800 -3.0070 1.1490 -10.4000 -1.489000
## 2000 meter buffer -1.20700 -2.5650 0.8169 -10.5900 -1.439000
## 1750 meter buffer -0.48880 -1.2580 1.0740 -4.9820 -1.055000
## 1500 meter buffer -0.47110 -0.6842 1.7070 -1.9960 -0.400400
## 500 meter buffer 0.85410 1.0740 2.2480 1.8080 0.939700
## 1250 meter buffer -1.29400 -1.3860 0.4148 -0.4861 -0.881700
## 1000 meter buffer -0.01823 -0.4209 1.2830 -0.9088 -0.007037
## 750 meter buffer -0.47550 -0.4547 0.2436 -0.9045 -0.389600
## cnd(ppl) cnd(rds) cnd(ssm_lns) cnd(trl) cnd(wll) df logLik
## 4500 meter buffer 7.559 4.0810 3.1700 -1.43600 2.705 15 -288.959
## 4750 meter buffer 6.605 3.8300 2.7500 -2.20800 4.448 15 -292.710
## 5000 meter buffer 5.895 2.9110 1.8420 -2.91200 3.735 15 -296.755
## 4250 meter buffer 4.318 2.5190 1.6940 -2.48800 4.988 15 -299.573
## 3750 meter buffer 4.796 4.2610 2.8580 -0.61380 6.808 15 -299.732
## 3500 meter buffer 3.627 3.6870 2.3370 -0.91640 7.028 15 -301.429
## 4000 meter buffer 3.917 3.7270 2.7440 -1.17400 5.774 15 -303.248
## 3000 meter buffer 4.395 2.1730 1.5150 -0.65670 5.858 15 -309.328
## 2750 meter buffer 4.360 2.5050 1.3730 -0.56820 6.085 15 -310.615
## 3250 meter buffer 4.167 2.5700 1.4920 -0.82350 4.721 15 -311.290
## 2500 meter buffer 4.748 3.1070 1.4460 0.41930 4.017 15 -314.347
## 2250 meter buffer 3.997 2.7960 0.8063 -0.02442 3.049 15 -317.267
## 2000 meter buffer 3.205 2.9030 0.7974 -0.18790 4.019 15 -320.493
## 1750 meter buffer 3.487 2.8000 0.8034 -0.11560 2.610 15 -325.585
## 1500 meter buffer 3.525 2.8130 0.9142 0.43940 2.328 15 -330.443
## 500 meter buffer 1.341 0.2240 -0.8317 -0.07879 1.632 15 -335.350
## 1250 meter buffer 3.134 2.1970 1.3760 0.41580 2.952 15 -337.221
## 1000 meter buffer 2.077 1.5330 0.2963 -0.39770 1.261 15 -340.337
## 750 meter buffer 1.779 0.7496 -0.2625 0.11670 1.559 15 -343.311
## AICc delta weight
## 4500 meter buffer 611.4 0.00 0.977
## 4750 meter buffer 619.0 7.50 0.023
## 5000 meter buffer 627.0 15.59 0.000
## 4250 meter buffer 632.7 21.23 0.000
## 3750 meter buffer 633.0 21.55 0.000
## 3500 meter buffer 636.4 24.94 0.000
## 4000 meter buffer 640.0 28.58 0.000
## 3000 meter buffer 652.2 40.74 0.000
## 2750 meter buffer 654.8 43.31 0.000
## 3250 meter buffer 656.1 44.66 0.000
## 2500 meter buffer 662.2 50.77 0.000
## 2250 meter buffer 668.1 56.62 0.000
## 2000 meter buffer 674.5 63.07 0.000
## 1750 meter buffer 684.7 73.25 0.000
## 1500 meter buffer 694.4 82.97 0.000
## 500 meter buffer 704.2 92.78 0.000
## 1250 meter buffer 708.0 96.52 0.000
## 1000 meter buffer 714.2 102.76 0.000
## 750 meter buffer 720.2 108.70 0.000
## Models ranked by AICc(x)
## Random terms (all models):
## cond(1 | array)for coyote top model appears to be 4500 m by quite a bit
Let’s get the model summary for this model
summary(coyote_mods_no250$`4500 meter buffer`)## Family: binomial ( logit )
## Formula:
## cbind(coyote, absent_coyote) ~ seismic_lines + pipeline + borrowpits +
## wellsites + roads + trails + lc_class20 + lc_class34 + lc_class50 +
## lc_class110 + lc_class210 + lc_class220 + lc_class230 + (1 | array)
## Data: .
##
## AIC BIC logLik deviance df.resid
## 607.9 653.3 -289.0 577.9 137
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## array (Intercept) 1.139 1.067
## Number of obs: 152, groups: array, 4
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.6021 1.9942 -0.302 0.762714
## seismic_lines 3.1702 1.6285 1.947 0.051570 .
## pipeline 7.5595 2.0938 3.610 0.000306 ***
## borrowpits -2.9490 6.7837 -0.435 0.663767
## wellsites 2.7046 3.4670 0.780 0.435331
## roads 4.0807 1.6767 2.434 0.014944 *
## trails -1.4364 1.6043 -0.895 0.370593
## lc_class20 -1.2981 2.2657 -0.573 0.566686
## lc_class34 -80.3309 26.9658 -2.979 0.002892 **
## lc_class50 -5.6739 1.9256 -2.947 0.003213 **
## lc_class110 -4.1190 2.8545 -1.443 0.149031
## lc_class210 -4.5776 1.9520 -2.345 0.019023 *
## lc_class220 -4.2835 2.6021 -1.646 0.099722 .
## lc_class230 4.2893 2.4831 1.727 0.084092 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1There is one lc class with a very high estimate and SE which seems a bit sus to me
Fisher
fisher_mods_no250 <- final_df %>%
# remove 350 meter buffer width
purrr::discard_at('250 meter buffer') %>%
# use purrr map to make global models for all other buffer sizes
purrr::map(
~.x %>%
glmmTMB::glmmTMB(cbind(fisher, absent_fisher) ~
seismic_lines +
pipeline +
borrowpits +
wellsites +
roads +
trails +
lc_class20 +
lc_class34 +
lc_class50 +
lc_class110 +
lc_class210 +
lc_class220 +
lc_class230 +
(1|array),
data = .,
family = 'binomial'))
# model selection
model.sel(fisher_mods_no250)## Model selection table
## cnd((Int)) dsp((Int)) cnd(brr) cnd(lc_c11) cnd(lc_c20)
## 2000 meter buffer 0.4777 + -32.590 -3.567 1.17400
## 2500 meter buffer 1.4910 + -38.920 -5.688 0.14320
## 3500 meter buffer 4.1770 + -30.180 -7.833 3.62600
## 2250 meter buffer 1.5800 + -28.730 -4.840 -0.22190
## 4250 meter buffer 3.9620 + -38.380 -5.465 -0.36770
## 3750 meter buffer 4.7270 + -33.610 -8.489 2.46100
## 2750 meter buffer 2.7180 + -28.750 -6.798 0.24470
## 1750 meter buffer -0.7824 + -25.720 -2.397 -0.87340
## 3000 meter buffer 2.8350 + -28.970 -5.665 2.19200
## 4500 meter buffer 3.6100 + -31.710 -3.062 1.49500
## 4750 meter buffer 4.4150 + -29.060 -6.770 2.15000
## 1500 meter buffer -1.0740 + -18.930 -2.022 -3.83300
## 3250 meter buffer 2.9140 + -25.880 -5.864 2.04000
## 4000 meter buffer 2.6360 + -34.130 -7.247 -0.04431
## 5000 meter buffer 6.3080 + -26.490 -11.850 0.37870
## 1000 meter buffer -1.2700 + -8.959 -2.474 -4.78300
## 750 meter buffer -1.3550 + -7.810 -1.681 -3.89700
## 1250 meter buffer -0.9845 + -12.050 -1.403 -3.05400
## 500 meter buffer -0.7114 + -4.084 -1.684 -3.46900
## cnd(lc_c21) cnd(lc_c22) cnd(lc_c23) cnd(lc_c34) cnd(lc_c50)
## 2000 meter buffer -2.6660 -0.4451 -1.3420 23.700 -2.744
## 2500 meter buffer -2.8280 -0.3343 -3.3040 34.530 -3.531
## 3500 meter buffer -4.6740 -3.2840 -3.9300 67.290 -4.821
## 2250 meter buffer -3.6800 -2.2220 -2.5310 21.600 -3.329
## 4250 meter buffer -7.4950 -3.3040 -2.0010 110.500 -4.232
## 3750 meter buffer -7.1490 -5.2890 -4.0780 55.500 -5.652
## 2750 meter buffer -3.8320 -2.8630 -3.6330 27.120 -4.572
## 1750 meter buffer -1.8020 0.2226 -2.2610 12.950 -2.622
## 3000 meter buffer -4.5900 -3.8870 -1.3630 36.910 -4.484
## 4500 meter buffer -9.0890 -4.6380 1.5650 101.500 -3.802
## 4750 meter buffer -11.0600 -5.1810 -0.7792 64.940 -6.207
## 1500 meter buffer -1.1410 0.9684 -3.1510 3.169 -2.077
## 3250 meter buffer -4.1410 -3.7680 -2.1750 34.060 -5.068
## 4000 meter buffer -6.4660 -1.8640 -4.5890 59.370 -4.925
## 5000 meter buffer -12.7400 -5.6570 -4.5260 49.940 -9.017
## 1000 meter buffer -2.3730 -1.4410 -1.2220 4.129 -1.552
## 750 meter buffer -0.9835 -0.9428 -1.4790 2.494 -1.228
## 1250 meter buffer -1.9370 -0.4350 -2.2990 2.347 -2.050
## 500 meter buffer -1.9700 -1.5380 -4.1550 -1.827 -2.673
## cnd(ppl) cnd(rds) cnd(ssm_lns) cnd(trl) cnd(wll) df logLik
## 2000 meter buffer 2.8820 -4.671 -2.85100 -2.34900 2.2470 15 -161.732
## 2500 meter buffer 2.0750 -4.868 -3.33900 -1.96200 3.1760 15 -162.116
## 3500 meter buffer 2.7680 -6.306 -5.11800 -3.91800 -2.4840 15 -164.414
## 2250 meter buffer 2.8520 -4.143 -2.77800 -2.12500 1.6660 15 -165.237
## 4250 meter buffer 1.9700 -6.669 -5.00400 -3.28200 -0.1988 15 -165.660
## 3750 meter buffer 1.4640 -4.996 -4.32100 -3.45900 1.9980 15 -166.455
## 2750 meter buffer 2.0280 -4.401 -3.50500 -2.13200 2.8010 15 -166.623
## 1750 meter buffer 0.9861 -3.050 -1.39000 -1.08900 5.4670 15 -166.963
## 3000 meter buffer 2.9530 -4.442 -4.11500 -2.67900 -1.3930 15 -167.753
## 4500 meter buffer 3.0610 -5.967 -5.14700 -3.67300 -5.3320 15 -168.688
## 4750 meter buffer 6.0540 -4.591 -3.59200 -2.89100 -5.9440 15 -168.903
## 1500 meter buffer 1.1980 -2.252 -1.49300 0.85810 1.2630 15 -168.993
## 3250 meter buffer 3.0550 -4.403 -3.93900 -2.49400 -0.9106 15 -169.349
## 4000 meter buffer 2.9760 -3.921 -2.44600 -1.43600 0.3342 15 -169.486
## 5000 meter buffer 5.7250 -4.580 -2.73800 -3.17000 -0.7984 15 -169.656
## 1000 meter buffer 0.8024 -1.702 -0.12900 0.04776 3.2900 15 -171.386
## 750 meter buffer 0.3745 -2.168 -0.76870 -0.20040 1.8520 15 -172.179
## 1250 meter buffer 0.9903 -2.099 -1.03800 0.20690 1.7370 15 -173.439
## 500 meter buffer 0.3115 -1.448 -0.03289 -0.23560 -0.1009 15 -175.153
## AICc delta weight
## 2000 meter buffer 357.0 0.00 0.547
## 2500 meter buffer 357.8 0.77 0.373
## 3500 meter buffer 362.4 5.36 0.037
## 2250 meter buffer 364.0 7.01 0.016
## 4250 meter buffer 364.8 7.86 0.011
## 3750 meter buffer 366.4 9.45 0.005
## 2750 meter buffer 366.8 9.78 0.004
## 1750 meter buffer 367.5 10.46 0.003
## 3000 meter buffer 369.0 12.04 0.001
## 4500 meter buffer 370.9 13.91 0.001
## 4750 meter buffer 371.3 14.34 0.000
## 1500 meter buffer 371.5 14.52 0.000
## 3250 meter buffer 372.2 15.23 0.000
## 4000 meter buffer 372.5 15.51 0.000
## 5000 meter buffer 372.8 15.85 0.000
## 1000 meter buffer 376.3 19.31 0.000
## 750 meter buffer 377.9 20.89 0.000
## 1250 meter buffer 380.4 23.41 0.000
## 500 meter buffer 383.8 26.84 0.000
## Models ranked by AICc(x)
## Random terms (all models):
## cond(1 | array)For fisher top model is 2000 meter
Let’s print the summary for this model
summary(fisher_mods_no250$`2000 meter buffer`)## Family: binomial ( logit )
## Formula:
## cbind(fisher, absent_fisher) ~ seismic_lines + pipeline + borrowpits +
## wellsites + roads + trails + lc_class20 + lc_class34 + lc_class50 +
## lc_class110 + lc_class210 + lc_class220 + lc_class230 + (1 | array)
## Data: .
##
## AIC BIC logLik deviance df.resid
## 353.5 398.8 -161.7 323.5 137
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## array (Intercept) 1.203 1.097
## Number of obs: 152, groups: array, 4
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.4777 1.6407 0.291 0.77093
## seismic_lines -2.8511 1.2653 -2.253 0.02424 *
## pipeline 2.8825 1.5424 1.869 0.06165 .
## borrowpits -32.5901 10.4890 -3.107 0.00189 **
## wellsites 2.2466 2.4280 0.925 0.35480
## roads -4.6707 1.5463 -3.021 0.00252 **
## trails -2.3490 1.3271 -1.770 0.07672 .
## lc_class20 1.1744 2.2860 0.514 0.60745
## lc_class34 23.7043 7.4238 3.193 0.00141 **
## lc_class50 -2.7441 1.7393 -1.578 0.11464
## lc_class110 -3.5668 2.4259 -1.470 0.14148
## lc_class210 -2.6659 1.8236 -1.462 0.14376
## lc_class220 -0.4451 2.3989 -0.186 0.85282
## lc_class230 -1.3420 2.1412 -0.627 0.53084
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Again lc_class34 has a very high standard error, we may not have enough data in this landcover class to use in the final analysis
Grey wolf
wolf_mods_no250 <- final_df %>%
# remove 350 meter buffer width
purrr::discard_at('250 meter buffer') %>%
# use purrr map to make global models for all other buffer sizes
purrr::map(
~.x %>%
glmmTMB::glmmTMB(cbind(grey_wolf, absent_grey_wolf) ~
seismic_lines +
pipeline +
borrowpits +
wellsites +
roads +
trails +
lc_class20 +
lc_class34 +
lc_class50 +
lc_class110 +
lc_class210 +
lc_class220 +
lc_class230 +
(1|array),
data = .,
family = 'binomial'))
# model selection
model.sel(wolf_mods_no250)## Model selection table
## cnd((Int)) dsp((Int)) cnd(brr) cnd(lc_c11) cnd(lc_c20)
## 4500 meter buffer -2.3760 + 7.4700 4.2190 3.48700
## 3750 meter buffer -1.1550 + 12.7600 1.7890 0.32770
## 4250 meter buffer -1.7210 + 6.0970 4.5270 5.16300
## 4000 meter buffer -2.3240 + 7.9510 2.8390 3.27200
## 3250 meter buffer -0.9283 + 17.7000 0.1020 -3.25800
## 4750 meter buffer -2.9270 + -0.8079 0.9542 4.28100
## 3500 meter buffer -0.4318 + 16.8700 0.9095 -2.00300
## 5000 meter buffer -2.4590 + -4.9270 1.5210 4.01300
## 3000 meter buffer -1.0830 + 13.7200 -1.8640 -3.14900
## 2750 meter buffer -0.2885 + 28.5500 -3.1310 -4.33000
## 2500 meter buffer -0.5582 + 27.9000 -4.7200 -2.88200
## 2250 meter buffer -1.5020 + 24.8200 -3.9260 -2.28900
## 1750 meter buffer -4.0630 + 9.2620 -1.7810 -0.60340
## 2000 meter buffer -3.4600 + 11.9100 -2.0070 -0.38070
## 1500 meter buffer -4.5960 + 2.6420 -2.2050 0.06206
## 1000 meter buffer -4.4760 + 13.5400 -1.0460 -0.46530
## 1250 meter buffer -3.8390 + 1.2320 -4.4190 -2.60800
## 500 meter buffer -3.2510 + -8.6130 0.8618 -1.14500
## 750 meter buffer -3.3800 + 9.9040 -2.5930 -1.06700
## cnd(lc_c21) cnd(lc_c22) cnd(lc_c23) cnd(lc_c34) cnd(lc_c50)
## 4500 meter buffer -2.3410 -4.0690 5.2810 -38.0100 -0.06745
## 3750 meter buffer -4.9110 -4.7840 2.3690 -3.8560 -1.02500
## 4250 meter buffer -1.3560 -4.4950 5.0770 9.2130 -1.21900
## 4000 meter buffer -1.6870 -3.7840 3.9670 10.2700 -1.05900
## 3250 meter buffer -5.5170 -4.1020 -0.5360 -32.8700 -1.73200
## 4750 meter buffer -3.8870 -2.9150 5.0500 -43.1100 -0.68490
## 3500 meter buffer -5.1760 -5.5610 0.5805 -39.7300 -1.21000
## 5000 meter buffer -4.1420 -4.4310 4.9280 -16.5800 -1.65800
## 3000 meter buffer -4.9150 -2.5200 -2.4300 -38.5500 -3.48300
## 2750 meter buffer -5.1130 -4.3840 -3.6100 -39.7300 -3.97800
## 2500 meter buffer -4.4590 -3.3030 -3.7270 -41.1500 -4.02500
## 2250 meter buffer -3.4720 -2.3400 -3.2370 -33.2200 -3.48000
## 1750 meter buffer -2.4700 1.4140 -1.2070 -16.0800 -2.08800
## 2000 meter buffer -1.8580 0.8223 -1.5850 -22.6500 -2.89700
## 1500 meter buffer -1.5350 0.9774 -0.7343 -8.5490 -2.06900
## 1000 meter buffer 0.4203 -0.5070 0.9747 7.9030 1.66400
## 1250 meter buffer -2.0490 -0.3537 -2.8630 -2.9000 -1.87800
## 500 meter buffer -0.2757 -0.8958 -1.1020 0.7220 -1.24600
## 750 meter buffer -1.6750 -2.0430 -1.2560 0.9455 -1.86500
## cnd(ppl) cnd(rds) cnd(ssm_lns) cnd(trl) cnd(wll) df logLik
## 4500 meter buffer -14.4600 0.8636 -0.51310 -3.620000 13.5700 15 -181.075
## 3750 meter buffer -11.5500 0.2491 0.05153 -2.376000 11.1700 15 -181.192
## 4250 meter buffer -13.9700 -0.5907 -1.49700 -4.575000 12.2900 15 -181.309
## 4000 meter buffer -12.5400 0.5537 -0.46370 -2.282000 10.7300 15 -181.769
## 3250 meter buffer -9.4940 1.0920 1.43100 0.006818 7.3050 15 -183.056
## 4750 meter buffer -12.2800 2.3290 0.65270 -2.341000 15.4100 15 -183.836
## 3500 meter buffer -9.3540 0.6707 0.44110 -1.106000 7.2190 15 -184.565
## 5000 meter buffer -11.9200 2.0740 0.62640 -2.012000 15.0300 15 -185.288
## 3000 meter buffer -7.2760 1.5810 2.07800 1.012000 8.7270 15 -186.071
## 2750 meter buffer -6.9390 1.4380 2.21200 1.648000 6.6220 15 -186.132
## 2500 meter buffer -6.1300 1.7950 2.10800 2.240000 5.6230 15 -188.075
## 2250 meter buffer -4.8320 2.3570 2.28300 2.694000 3.8210 15 -190.528
## 1750 meter buffer -3.2660 3.6090 3.02700 2.463000 4.7360 15 -191.657
## 2000 meter buffer -3.6110 3.1190 2.55000 2.396000 4.5130 15 -191.954
## 1500 meter buffer -2.2150 4.2550 3.26900 3.050000 4.5690 15 -192.541
## 1000 meter buffer 0.9378 -0.0308 1.57800 2.127000 -1.9680 15 -193.218
## 1250 meter buffer -0.6217 3.8050 3.93400 3.359000 3.1410 15 -194.174
## 500 meter buffer 0.8645 1.4090 1.17900 1.436000 -1.0810 15 -199.340
## 750 meter buffer 1.9910 2.0850 2.72600 2.548000 0.4032 15 -199.871
## AICc delta weight
## 4500 meter buffer 395.7 0.00 0.291
## 3750 meter buffer 395.9 0.23 0.259
## 4250 meter buffer 396.1 0.47 0.230
## 4000 meter buffer 397.1 1.39 0.145
## 3250 meter buffer 399.6 3.96 0.040
## 4750 meter buffer 401.2 5.52 0.018
## 3500 meter buffer 402.7 6.98 0.009
## 5000 meter buffer 404.1 8.43 0.004
## 3000 meter buffer 405.7 9.99 0.002
## 2750 meter buffer 405.8 10.11 0.002
## 2500 meter buffer 409.7 14.00 0.000
## 2250 meter buffer 414.6 18.90 0.000
## 1750 meter buffer 416.8 21.16 0.000
## 2000 meter buffer 417.4 21.76 0.000
## 1500 meter buffer 418.6 22.93 0.000
## 1000 meter buffer 420.0 24.29 0.000
## 1250 meter buffer 421.9 26.20 0.000
## 500 meter buffer 432.2 36.53 0.000
## 750 meter buffer 433.3 37.59 0.000
## Models ranked by AICc(x)
## Random terms (all models):
## cond(1 | array)For grey wolf top model is 4500 m buffer
Let’s get the model summary for this model
summary(wolf_mods_no250$`4500 meter buffer`)## Family: binomial ( logit )
## Formula:
## cbind(grey_wolf, absent_grey_wolf) ~ seismic_lines + pipeline +
## borrowpits + wellsites + roads + trails + lc_class20 + lc_class34 +
## lc_class50 + lc_class110 + lc_class210 + lc_class220 + lc_class230 +
## (1 | array)
## Data: .
##
## AIC BIC logLik deviance df.resid
## 392.2 437.5 -181.1 362.2 137
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## array (Intercept) 3.88e-09 6.229e-05
## Number of obs: 152, groups: array, 4
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -2.37568 2.27310 -1.045 0.295964
## seismic_lines -0.51313 1.61209 -0.318 0.750257
## pipeline -14.45950 3.67826 -3.931 8.46e-05 ***
## borrowpits 7.46992 9.41055 0.794 0.427322
## wellsites 13.56596 3.67755 3.689 0.000225 ***
## roads 0.86359 1.88444 0.458 0.646757
## trails -3.62005 2.08054 -1.740 0.081867 .
## lc_class20 3.48719 3.38830 1.029 0.303392
## lc_class34 -38.00581 37.06030 -1.026 0.305121
## lc_class50 -0.06745 2.55025 -0.026 0.978901
## lc_class110 4.21886 3.56984 1.182 0.237282
## lc_class210 -2.34074 2.99296 -0.782 0.434166
## lc_class220 -4.06945 3.58435 -1.135 0.256233
## lc_class230 5.28117 3.03467 1.740 0.081810 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1lc_class34 still presenting some issues, interesting that seismic lines weren’t significant and have a negative estimate
Lynx
lynx_mods_no250 <- final_df %>%
# remove 350 meter buffer width
purrr::discard_at('250 meter buffer') %>%
# use purrr map to make global models for all other buffer sizes
purrr::map(
~.x %>%
glmmTMB::glmmTMB(cbind(lynx, absent_lynx) ~
seismic_lines +
pipeline +
borrowpits +
wellsites +
roads +
trails +
lc_class20 +
lc_class34 +
lc_class50 +
lc_class110 +
lc_class210 +
lc_class220 +
lc_class230 +
(1|array),
data = .,
family = 'binomial'))
# model selection
model.sel(lynx_mods_no250)## Model selection table
## cnd((Int)) dsp((Int)) cnd(brr) cnd(lc_c11) cnd(lc_c20)
## 1000 meter buffer -2.2210 + 5.399 -1.95100 0.8875
## 2750 meter buffer -2.0820 + -14.450 -2.65100 0.6012
## 1250 meter buffer -3.4560 + -2.301 -1.40500 0.8151
## 2500 meter buffer -4.6450 + -10.490 -0.82550 1.4370
## 3000 meter buffer -1.6070 + -7.243 -3.53000 -1.6480
## 3250 meter buffer -0.9863 + -8.403 -3.23000 0.4247
## 4750 meter buffer -0.9497 + -8.377 -8.63200 -2.2350
## 750 meter buffer -2.3060 + 4.079 -1.38700 0.1944
## 5000 meter buffer -0.5495 + -7.298 -6.77600 -2.8140
## 1500 meter buffer -2.9580 + 1.707 -2.85200 0.2368
## 4000 meter buffer -0.1394 + -6.980 -6.92800 -1.1010
## 3750 meter buffer -0.6332 + -5.122 -6.14900 -1.9300
## 3500 meter buffer -0.3969 + -3.435 -4.72400 -0.8918
## 4500 meter buffer -0.7145 + -2.491 -6.03000 -2.2220
## 2000 meter buffer -3.2160 + 1.999 -1.60400 -0.4112
## 4250 meter buffer -0.3378 + -3.711 -5.32300 -1.9140
## 500 meter buffer -3.3190 + -6.902 -0.80020 0.2681
## 2250 meter buffer -4.2780 + -3.953 0.03378 2.1620
## 1750 meter buffer -2.8810 + 0.271 -2.60000 0.1709
## cnd(lc_c21) cnd(lc_c22) cnd(lc_c23) cnd(lc_c34) cnd(lc_c50)
## 1000 meter buffer 1.390000 -2.7410 1.2020 -7.307 -1.5690
## 2750 meter buffer 0.201500 -0.4946 -1.7070 -32.490 -4.7990
## 1250 meter buffer 1.800000 -1.2450 1.4820 -8.094 -0.6688
## 2500 meter buffer 0.661300 0.7638 0.7968 -29.420 -2.8820
## 3000 meter buffer -0.198300 -0.9439 -3.1130 -42.730 -4.4660
## 3250 meter buffer 0.004042 -1.9630 -2.0430 -30.450 -4.1350
## 4750 meter buffer 1.287000 -1.1540 -6.2360 -51.700 -4.0210
## 750 meter buffer 1.156000 -1.8820 0.8398 -3.056 -0.9768
## 5000 meter buffer 2.050000 -1.9720 -5.2110 -12.270 -4.3590
## 1500 meter buffer 0.168300 -2.1120 0.6480 -11.400 -1.5960
## 4000 meter buffer 0.234100 -2.9940 -5.3570 -39.620 -5.0820
## 3750 meter buffer -0.264200 -2.9340 -4.2780 -44.570 -4.5390
## 3500 meter buffer -0.139500 -3.0250 -3.4490 -32.640 -4.3190
## 4500 meter buffer 0.969500 -2.6250 -5.3680 -58.530 -3.2620
## 2000 meter buffer -0.021910 -1.8120 -0.1756 -24.210 -2.6230
## 4250 meter buffer 0.681700 -3.2270 -5.0030 -46.610 -4.0660
## 500 meter buffer 1.112000 0.3428 2.8550 1.288 0.6173
## 2250 meter buffer 0.454400 -0.4712 1.3290 -20.490 -1.5680
## 1750 meter buffer 0.770000 -1.8290 -0.4908 -12.120 -2.1160
## cnd(ppl) cnd(rds) cnd(ssm_lns) cnd(trl) cnd(wll) df logLik
## 1000 meter buffer 2.4440 3.281 0.4976 1.20000 -1.7090 15 -275.660
## 2750 meter buffer 2.3410 3.351 2.1430 -0.02248 9.8790 15 -285.193
## 1250 meter buffer 2.9570 3.928 1.1710 2.02300 -0.2275 15 -288.524
## 2500 meter buffer 2.1260 5.256 3.3280 2.21200 7.2750 15 -288.859
## 3000 meter buffer 2.2080 3.087 2.5960 0.51900 9.1360 15 -289.307
## 3250 meter buffer 0.5244 1.960 1.5450 -1.10100 9.9450 15 -290.161
## 4750 meter buffer -0.7878 2.946 3.0040 0.73550 14.4400 15 -290.493
## 750 meter buffer 1.8540 2.377 0.4768 0.94980 -1.1550 15 -291.458
## 5000 meter buffer -2.6820 1.680 2.1100 -0.13930 14.8300 15 -291.620
## 1500 meter buffer 2.2640 4.107 1.9640 1.98300 2.8990 15 -292.060
## 4000 meter buffer 1.5200 2.162 2.6190 -0.18940 10.1100 15 -292.111
## 3750 meter buffer 1.1240 2.788 2.8460 0.22840 10.5000 15 -292.247
## 3500 meter buffer 0.6192 1.829 1.9490 -0.82650 9.8870 15 -292.580
## 4500 meter buffer 0.5021 2.481 2.4800 0.28500 10.4200 15 -294.536
## 2000 meter buffer 0.7541 4.724 2.8600 2.44200 4.4220 15 -294.941
## 4250 meter buffer 1.5860 2.010 2.4000 -0.11440 8.0520 15 -295.023
## 500 meter buffer 0.6402 1.696 0.4167 0.21630 2.3700 15 -295.864
## 2250 meter buffer 2.3120 4.562 2.3530 2.11700 4.4910 15 -296.883
## 1750 meter buffer 1.5140 3.823 2.1730 1.91100 3.6580 15 -297.309
## AICc delta weight
## 1000 meter buffer 584.9 0.00 1
## 2750 meter buffer 603.9 19.06 0
## 1250 meter buffer 610.6 25.73 0
## 2500 meter buffer 611.2 26.40 0
## 3000 meter buffer 612.1 27.29 0
## 3250 meter buffer 613.9 29.00 0
## 4750 meter buffer 614.5 29.66 0
## 750 meter buffer 616.4 31.59 0
## 5000 meter buffer 616.8 31.92 0
## 1500 meter buffer 617.6 32.80 0
## 4000 meter buffer 617.8 32.90 0
## 3750 meter buffer 618.0 33.17 0
## 3500 meter buffer 618.7 33.84 0
## 4500 meter buffer 622.6 37.75 0
## 2000 meter buffer 623.4 38.56 0
## 4250 meter buffer 623.6 38.73 0
## 500 meter buffer 625.3 40.41 0
## 2250 meter buffer 627.3 42.45 0
## 1750 meter buffer 628.1 43.30 0
## Models ranked by AICc(x)
## Random terms (all models):
## cond(1 | array)For lynx the top model is the 1000 m buffer
Let’s get the model summary
summary(lynx_mods_no250$`1000 meter buffer`)## Family: binomial ( logit )
## Formula:
## cbind(lynx, absent_lynx) ~ seismic_lines + pipeline + borrowpits +
## wellsites + roads + trails + lc_class20 + lc_class34 + lc_class50 +
## lc_class110 + lc_class210 + lc_class220 + lc_class230 + (1 | array)
## Data: .
##
## AIC BIC logLik deviance df.resid
## 581.3 626.7 -275.7 551.3 137
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## array (Intercept) 0.005282 0.07268
## Number of obs: 152, groups: array, 4
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -2.2213 1.7785 -1.249 0.21167
## seismic_lines 0.4976 0.6728 0.740 0.45959
## pipeline 2.4438 0.9253 2.641 0.00826 **
## borrowpits 5.3992 6.6173 0.816 0.41455
## wellsites -1.7095 1.5583 -1.097 0.27265
## roads 3.2815 0.8106 4.048 5.16e-05 ***
## trails 1.2002 0.7234 1.659 0.09711 .
## lc_class20 0.8875 1.8019 0.493 0.62235
## lc_class34 -7.3069 5.3934 -1.355 0.17549
## lc_class50 -1.5694 1.9859 -0.790 0.42937
## lc_class110 -1.9508 1.9206 -1.016 0.30975
## lc_class210 1.3900 1.7371 0.800 0.42360
## lc_class220 -2.7407 2.1055 -1.302 0.19303
## lc_class230 1.2019 2.0790 0.578 0.56319
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Moose
moose_mods_no250 <- final_df %>%
# remove 350 meter buffer width
purrr::discard_at('250 meter buffer') %>%
# use purrr map to make global models for all other buffer sizes
purrr::map(
~.x %>%
glmmTMB::glmmTMB(cbind(moose, absent_moose) ~
seismic_lines +
pipeline +
borrowpits +
wellsites +
roads +
trails +
lc_class20 +
lc_class34 +
lc_class50 +
lc_class110 +
lc_class210 +
lc_class220 +
lc_class230 +
(1|array),
data = .,
family = 'binomial'))
# model selection
model.sel(moose_mods_no250)## Model selection table
## cnd((Int)) dsp((Int)) cnd(brr) cnd(lc_c11) cnd(lc_c20)
## 3750 meter buffer -9.205 + -11.8000 9.427 8.205
## 4000 meter buffer -9.283 + -14.4400 10.120 9.262
## 3250 meter buffer -7.326 + -0.6449 8.858 6.573
## 3000 meter buffer -5.532 + 1.4540 5.762 3.001
## 4250 meter buffer -9.420 + -15.3200 9.725 10.260
## 5000 meter buffer -9.865 + -14.2400 7.416 8.885
## 4500 meter buffer -9.886 + -16.8800 8.457 9.441
## 4750 meter buffer -10.090 + -15.7600 9.437 9.390
## 3500 meter buffer -8.396 + -6.9660 8.888 7.376
## 2750 meter buffer -5.325 + -0.4263 5.463 3.785
## 2250 meter buffer -4.683 + 2.5550 5.271 3.260
## 2000 meter buffer -5.559 + -4.1900 5.645 4.327
## 2500 meter buffer -4.989 + 2.2070 5.165 3.506
## 1750 meter buffer -4.634 + -4.4430 3.992 3.557
## 1250 meter buffer -3.377 + -6.5860 2.341 2.767
## 750 meter buffer -2.671 + -3.7510 2.683 2.479
## 1500 meter buffer -3.983 + -3.3340 3.989 3.356
## 500 meter buffer -2.884 + -9.3930 3.598 2.436
## 1000 meter buffer -2.632 + -1.5350 2.143 1.988
## cnd(lc_c21) cnd(lc_c22) cnd(lc_c23) cnd(lc_c34) cnd(lc_c50)
## 3750 meter buffer 8.258 9.361 7.209 -8.7110 7.035
## 4000 meter buffer 9.171 9.396 7.292 13.8500 6.964
## 3250 meter buffer 7.480 5.983 5.986 -13.9900 6.544
## 3000 meter buffer 6.201 5.076 3.220 -18.0600 4.122
## 4250 meter buffer 9.075 10.180 7.117 8.6850 6.945
## 5000 meter buffer 9.099 11.000 5.204 -35.2300 7.330
## 4500 meter buffer 8.458 10.870 6.332 -14.2200 7.256
## 4750 meter buffer 9.182 10.490 6.977 -16.1700 7.824
## 3500 meter buffer 7.840 7.591 6.965 -6.6260 6.840
## 2750 meter buffer 6.430 5.300 3.385 -6.8180 4.115
## 2250 meter buffer 5.627 4.575 3.495 -9.1220 3.350
## 2000 meter buffer 6.072 5.814 4.648 -4.8470 3.932
## 2500 meter buffer 5.667 4.218 3.464 -8.1650 3.508
## 1750 meter buffer 5.873 4.918 4.304 2.0300 4.677
## 1250 meter buffer 4.482 4.531 2.309 -1.3720 3.378
## 750 meter buffer 2.211 3.210 2.779 -5.4720 2.815
## 1500 meter buffer 4.602 4.499 3.761 -3.7120 3.564
## 500 meter buffer 2.934 3.223 2.892 0.8086 3.270
## 1000 meter buffer 2.757 2.122 3.182 -0.4888 2.968
## cnd(ppl) cnd(rds) cnd(ssm_lns) cnd(trl) cnd(wll) df logLik
## 3750 meter buffer -6.414 1.42300 0.82500 0.2511 -2.9620 15 -288.239
## 4000 meter buffer -6.958 0.71860 0.48330 -0.1971 -2.7700 15 -289.516
## 3250 meter buffer -7.524 0.62500 0.08454 -0.5983 -1.9880 15 -291.215
## 3000 meter buffer -7.958 0.19990 0.44840 -0.6199 1.9760 15 -291.674
## 4250 meter buffer -5.966 0.69870 0.50480 -0.3154 -3.0590 15 -291.707
## 5000 meter buffer -4.758 1.86500 1.95400 1.0270 -2.0550 15 -292.065
## 4500 meter buffer -4.990 1.81200 1.53400 0.6604 -2.4790 15 -292.194
## 4750 meter buffer -5.872 1.70900 1.28300 0.5059 -1.9210 15 -292.959
## 3500 meter buffer -7.586 1.18400 0.57060 -0.2578 -1.2150 15 -292.978
## 2750 meter buffer -7.448 -0.28040 -0.18940 -1.1490 1.6850 15 -296.355
## 2250 meter buffer -6.493 -0.49530 -0.38530 -1.2870 0.4424 15 -298.679
## 2000 meter buffer -6.048 0.15560 -0.27670 -1.1850 -0.4839 15 -298.782
## 2500 meter buffer -6.713 0.04149 0.04839 -0.8332 0.4888 15 -298.965
## 1750 meter buffer -5.739 -0.77580 -1.04600 -1.3770 -2.2170 15 -299.799
## 1250 meter buffer -4.047 -1.06700 -1.17000 -1.5400 -3.2830 15 -302.520
## 750 meter buffer -3.300 0.16850 -1.19400 -2.1610 -1.1880 15 -302.567
## 1500 meter buffer -5.221 -0.72240 -1.09500 -1.5290 -1.8400 15 -303.054
## 500 meter buffer -2.982 -0.88640 -1.61300 -2.0210 -0.7721 15 -307.616
## 1000 meter buffer -3.031 -1.04500 -1.14000 -1.7570 -1.6820 15 -310.178
## AICc delta weight
## 3750 meter buffer 610.0 0.00 0.689
## 4000 meter buffer 612.6 2.55 0.192
## 3250 meter buffer 616.0 5.95 0.035
## 3000 meter buffer 616.9 6.87 0.022
## 4250 meter buffer 616.9 6.94 0.021
## 5000 meter buffer 617.7 7.65 0.015
## 4500 meter buffer 617.9 7.91 0.013
## 4750 meter buffer 619.4 9.44 0.006
## 3500 meter buffer 619.5 9.48 0.006
## 2750 meter buffer 626.2 16.23 0.000
## 2250 meter buffer 630.9 20.88 0.000
## 2000 meter buffer 631.1 21.09 0.000
## 2500 meter buffer 631.5 21.45 0.000
## 1750 meter buffer 633.1 23.12 0.000
## 1250 meter buffer 638.6 28.56 0.000
## 750 meter buffer 638.7 28.66 0.000
## 1500 meter buffer 639.6 29.63 0.000
## 500 meter buffer 648.8 38.75 0.000
## 1000 meter buffer 653.9 43.88 0.000
## Models ranked by AICc(x)
## Random terms (all models):
## cond(1 | array)For moose the top model is the 3750 m buffer
Let’s get the model summary
summary(moose_mods_no250$`3750 meter buffer`)## Family: binomial ( logit )
## Formula:
## cbind(moose, absent_moose) ~ seismic_lines + pipeline + borrowpits +
## wellsites + roads + trails + lc_class20 + lc_class34 + lc_class50 +
## lc_class110 + lc_class210 + lc_class220 + lc_class230 + (1 | array)
## Data: .
##
## AIC BIC logLik deviance df.resid
## 606.5 651.8 -288.2 576.5 137
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## array (Intercept) 2.321e-10 1.524e-05
## Number of obs: 152, groups: array, 4
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -9.2052 1.5989 -5.757 8.56e-09 ***
## seismic_lines 0.8250 0.8882 0.929 0.352960
## pipeline -6.4141 1.6781 -3.822 0.000132 ***
## borrowpits -11.8016 6.0566 -1.949 0.051351 .
## wellsites -2.9624 2.1366 -1.387 0.165593
## roads 1.4234 1.0848 1.312 0.189470
## trails 0.2511 1.1182 0.225 0.822354
## lc_class20 8.2051 2.1166 3.877 0.000106 ***
## lc_class34 -8.7109 14.3630 -0.606 0.544196
## lc_class50 7.0351 1.7421 4.038 5.38e-05 ***
## lc_class110 9.4270 2.0348 4.633 3.61e-06 ***
## lc_class210 8.2584 1.7866 4.622 3.79e-06 ***
## lc_class220 9.3615 2.1734 4.307 1.65e-05 ***
## lc_class230 7.2091 1.8109 3.981 6.86e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Red fox
fox_mods_no250 <- final_df %>%
# remove 350 meter buffer width
purrr::discard_at('250 meter buffer') %>%
# use purrr map to make global models for all other buffer sizes
purrr::map(
~.x %>%
glmmTMB::glmmTMB(cbind(red_fox, absent_red_fox) ~
seismic_lines +
pipeline +
borrowpits +
wellsites +
roads +
trails +
lc_class20 +
lc_class34 +
lc_class50 +
lc_class110 +
lc_class210 +
lc_class220 +
lc_class230 +
(1|array),
data = .,
family = 'binomial'))
# model selection
model.sel(fox_mods_no250)## Model selection table
## cnd((Int)) dsp((Int)) cnd(brr) cnd(lc_c11) cnd(lc_c20)
## 3750 meter buffer -23.690 + -119.700 24.54000 8.1590
## 4000 meter buffer -27.740 + -79.470 30.62000 7.5130
## 500 meter buffer -6.447 + 14.000 0.07924 1.1800
## 4750 meter buffer -33.510 + -104.300 28.70000 15.8700
## 4500 meter buffer -26.500 + -95.040 30.18000 6.5690
## 3500 meter buffer -18.740 + -85.750 19.81000 4.3940
## 4250 meter buffer -24.140 + -85.330 24.79000 9.8620
## 1000 meter buffer -9.402 + -10.500 3.56700 2.3310
## 3250 meter buffer -15.220 + -83.780 16.48000 4.3700
## 5000 meter buffer -27.550 + -101.300 20.31000 10.0200
## 750 meter buffer -5.662 + -7.859 0.40800 -1.4200
## 2250 meter buffer -5.211 + -62.400 -1.68000 1.4500
## 1250 meter buffer -7.659 + -7.811 0.87910 -1.7010
## 3000 meter buffer -10.390 + -91.500 8.56500 3.3110
## 2750 meter buffer -4.720 + -68.670 0.26310 0.5898
## 2500 meter buffer -5.085 + -52.980 0.54390 2.1990
## 2000 meter buffer -6.283 + -39.650 0.97650 2.6660
## 1500 meter buffer -8.468 + -6.456 4.25300 1.1330
## 1750 meter buffer -5.487 + -17.630 1.99600 1.4570
## cnd(lc_c21) cnd(lc_c22) cnd(lc_c23) cnd(lc_c34) cnd(lc_c50)
## 3750 meter buffer 10.5100 28.0700 16.600 110.9000 26.2400
## 4000 meter buffer 11.3200 29.2200 19.630 33.5900 27.8400
## 500 meter buffer 0.2117 -0.2819 -1.952 -3.2380 -1.9700
## 4750 meter buffer 12.5800 41.1700 22.280 187.4000 27.2100
## 4500 meter buffer 9.4540 25.5400 23.270 152.1000 24.1000
## 3500 meter buffer 9.9510 20.4000 14.760 56.4100 19.2600
## 4250 meter buffer 10.0700 30.1600 15.550 137.9000 24.2300
## 1000 meter buffer 6.7910 8.6490 2.478 -5.7110 1.2590
## 3250 meter buffer 8.5530 16.2400 13.760 41.4300 14.9200
## 5000 meter buffer 6.6050 35.4800 15.660 167.5000 22.1000
## 750 meter buffer 2.3640 2.0760 -1.931 -2.3510 -1.3010
## 2250 meter buffer 5.9970 4.9570 4.426 42.9600 7.2900
## 1250 meter buffer 5.3800 7.7410 -3.171 -7.1840 0.6127
## 3000 meter buffer 6.2900 9.6700 11.380 47.1100 12.0100
## 2750 meter buffer 4.9190 5.4770 3.450 50.8700 5.6880
## 2500 meter buffer 5.4580 5.4220 3.557 48.4900 5.2580
## 2000 meter buffer 5.8460 3.5100 4.655 26.1700 4.5180
## 1500 meter buffer 7.0080 7.3430 1.354 -0.8072 2.6440
## 1750 meter buffer 4.8150 1.7860 1.280 4.0420 1.6860
## cnd(ppl) cnd(rds) cnd(ssm_lns) cnd(trl) cnd(wll) df logLik
## 3750 meter buffer 3.09700 1.1080 -0.24620 4.4620 -19.4100 15 -128.145
## 4000 meter buffer 7.57700 5.6120 2.89800 7.1960 -28.7900 15 -128.278
## 500 meter buffer 5.02200 5.4990 3.21200 2.3900 -11.7200 15 -128.681
## 4750 meter buffer 0.03691 8.3610 5.63100 5.5740 -10.8900 15 -132.319
## 4500 meter buffer -2.93400 5.4840 3.20500 3.2140 -9.8480 15 -132.586
## 3500 meter buffer 6.77800 -0.1700 -0.83030 2.7110 -18.3400 15 -132.688
## 4250 meter buffer 2.37100 1.6930 0.03921 3.4580 -16.6100 15 -132.874
## 1000 meter buffer 1.47200 4.7540 0.17280 -0.5835 2.4390 15 -133.355
## 3250 meter buffer 6.66300 -1.5920 -2.27300 1.3510 -13.3800 15 -134.939
## 5000 meter buffer -2.80900 7.7850 6.27300 3.8360 -1.5170 15 -135.981
## 750 meter buffer 2.86900 3.2770 0.86510 0.5270 1.5320 15 -137.015
## 2250 meter buffer 8.16500 -4.4860 -4.94000 -1.8550 -8.0530 15 -138.154
## 1250 meter buffer -0.25240 3.3500 0.71720 2.3310 -0.5302 15 -138.835
## 3000 meter buffer 6.31500 -2.9290 -3.92700 0.3018 -8.9350 15 -139.954
## 2750 meter buffer 6.13300 -5.1950 -4.67000 -1.9830 -4.6950 15 -140.819
## 2500 meter buffer 4.89200 -4.7870 -4.53000 -2.4000 -2.8530 15 -142.155
## 2000 meter buffer 7.60400 -1.8660 -3.26900 0.3092 -7.7120 15 -142.797
## 1500 meter buffer 0.99900 1.1700 -1.19400 3.0700 -3.0510 15 -146.669
## 1750 meter buffer 1.01900 0.2511 -1.43300 0.5256 -0.4095 15 -153.047
## AICc delta weight
## 3750 meter buffer 289.8 0.00 0.398
## 4000 meter buffer 290.1 0.27 0.348
## 500 meter buffer 290.9 1.07 0.233
## 4750 meter buffer 298.2 8.35 0.006
## 4500 meter buffer 298.7 8.88 0.005
## 3500 meter buffer 298.9 9.08 0.004
## 4250 meter buffer 299.3 9.46 0.004
## 1000 meter buffer 300.2 10.42 0.002
## 3250 meter buffer 303.4 13.59 0.000
## 5000 meter buffer 305.5 15.67 0.000
## 750 meter buffer 307.6 17.74 0.000
## 2250 meter buffer 309.8 20.02 0.000
## 1250 meter buffer 311.2 21.38 0.000
## 3000 meter buffer 313.4 23.62 0.000
## 2750 meter buffer 315.2 25.35 0.000
## 2500 meter buffer 317.8 28.02 0.000
## 2000 meter buffer 319.1 29.30 0.000
## 1500 meter buffer 326.9 37.05 0.000
## 1750 meter buffer 339.6 49.80 0.000
## Models ranked by AICc(x)
## Random terms (all models):
## cond(1 | array)For red fox the top model is 3750 m buffer
Let’s get the model summary
summary(fox_mods_no250$`3750 meter buffer`)## Family: binomial ( logit )
## Formula:
## cbind(red_fox, absent_red_fox) ~ seismic_lines + pipeline + borrowpits +
## wellsites + roads + trails + lc_class20 + lc_class34 + lc_class50 +
## lc_class110 + lc_class210 + lc_class220 + lc_class230 + (1 | array)
## Data: .
##
## AIC BIC logLik deviance df.resid
## 286.3 331.6 -128.1 256.3 137
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## array (Intercept) 3.94 1.985
## Number of obs: 152, groups: array, 4
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -23.6900 5.9704 -3.968 7.25e-05 ***
## seismic_lines -0.2462 3.0527 -0.081 0.935711
## pipeline 3.0969 3.9456 0.785 0.432510
## borrowpits -119.6821 27.1704 -4.405 1.06e-05 ***
## wellsites -19.4084 7.8584 -2.470 0.013520 *
## roads 1.1082 3.5675 0.311 0.756080
## trails 4.4624 3.3651 1.326 0.184814
## lc_class20 8.1589 6.4236 1.270 0.204031
## lc_class34 110.9227 37.6596 2.945 0.003225 **
## lc_class50 26.2421 6.3973 4.102 4.09e-05 ***
## lc_class110 24.5357 6.6335 3.699 0.000217 ***
## lc_class210 10.5121 5.3219 1.975 0.048239 *
## lc_class220 28.0743 7.6920 3.650 0.000262 ***
## lc_class230 16.5970 6.5277 2.543 0.011005 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Gagh! Borrow pits does not have a reasonable estimate and SE
White-tailed deer
deer_mods_no250 <- final_df %>%
# remove 350 meter buffer width
purrr::discard_at('250 meter buffer') %>%
# use purrr map to make global models for all other buffer sizes
purrr::map(
~.x %>%
# have to include the `` around the white-tailed_deer or R won't recognize it as a variable because of the -
glmmTMB::glmmTMB(cbind(`white-tailed_deer`, `absent_white-tailed_deer`) ~
seismic_lines +
pipeline +
borrowpits +
wellsites +
roads +
trails +
lc_class20 +
lc_class34 +
lc_class50 +
lc_class110 +
lc_class210 +
lc_class220 +
lc_class230 +
(1|array),
data = .,
family = 'binomial'))
# model selection
model.sel(deer_mods_no250)## Model selection table
## cnd((Int)) dsp((Int)) cnd(brr) cnd(lc_c11) cnd(lc_c20)
## 3750 meter buffer 6.8600 + -16.910 -11.390 -6.093
## 4500 meter buffer 7.7950 + -18.770 -11.500 -6.895
## 3500 meter buffer 7.6200 + -9.968 -12.200 -8.638
## 5000 meter buffer 8.3290 + -21.840 -13.520 -7.094
## 4750 meter buffer 9.4510 + -18.680 -13.240 -7.075
## 4000 meter buffer 6.1040 + -22.400 -10.180 -5.315
## 4250 meter buffer 6.7000 + -20.320 -10.370 -4.932
## 2750 meter buffer 5.3670 + -14.960 -8.832 -7.014
## 3000 meter buffer 3.9600 + -18.760 -7.638 -5.160
## 2500 meter buffer 4.8710 + -14.020 -9.047 -7.365
## 3250 meter buffer 5.1720 + -10.120 -8.816 -6.733
## 2250 meter buffer 4.9470 + -9.646 -8.071 -6.468
## 2000 meter buffer 3.3110 + -16.950 -6.580 -5.965
## 1750 meter buffer 2.8700 + -16.870 -5.833 -4.865
## 1250 meter buffer 2.6220 + -14.500 -4.975 -3.611
## 1500 meter buffer 2.3030 + -15.890 -4.740 -4.111
## 500 meter buffer 0.9442 + -1.352 -1.340 -1.550
## 750 meter buffer 1.8810 + -3.276 -2.500 -1.737
## 1000 meter buffer 1.9330 + -11.230 -3.278 -2.742
## cnd(lc_c21) cnd(lc_c22) cnd(lc_c23) cnd(lc_c34) cnd(lc_c50)
## 3750 meter buffer -7.0520 -14.0300 -5.21700 -28.0300 -9.7310
## 4500 meter buffer -8.5460 -16.9000 -4.61600 -71.0900 -9.7670
## 3500 meter buffer -6.7340 -14.2300 -7.68400 -27.1100 -9.7160
## 5000 meter buffer -8.6500 -17.1000 -5.28800 -57.6700 -11.4900
## 4750 meter buffer -9.1870 -17.9500 -5.97200 -51.6100 -11.5000
## 4000 meter buffer -6.1680 -12.4100 -4.20800 -4.2870 -9.1920
## 4250 meter buffer -7.0870 -13.3000 -3.52700 -24.6900 -8.8920
## 2750 meter buffer -4.2370 -10.3100 -4.99900 4.0600 -7.1010
## 3000 meter buffer -4.7470 -8.8680 -2.60300 -4.3240 -6.4930
## 2500 meter buffer -3.5460 -8.7360 -6.30100 4.8390 -7.4870
## 3250 meter buffer -5.7660 -10.9700 -4.38500 -16.8600 -7.0020
## 2250 meter buffer -3.7260 -8.3040 -5.59700 -0.9146 -6.6250
## 2000 meter buffer -2.9360 -5.6380 -3.42400 6.1230 -5.1200
## 1750 meter buffer -2.7070 -4.5180 -3.89800 6.5230 -4.7930
## 1250 meter buffer -2.1320 -3.8580 -4.38400 6.0270 -3.2730
## 1500 meter buffer -2.7060 -3.8610 -3.85100 0.4099 -4.3870
## 500 meter buffer -0.5153 -0.4099 0.07816 0.4883 -0.2481
## 750 meter buffer -1.8120 -1.1940 -2.01500 -2.2460 -2.3280
## 1000 meter buffer -1.8140 -2.7900 -1.87900 2.7100 -2.2840
## cnd(ppl) cnd(rds) cnd(ssm_lns) cnd(trl) cnd(wll) df logLik
## 3750 meter buffer 2.5130 2.7930 2.1670 0.1252 6.99600 15 -350.163
## 4500 meter buffer 5.5380 3.4430 2.4770 0.4691 1.64300 15 -351.755
## 3500 meter buffer 1.4700 2.2170 2.1800 0.2082 7.59500 15 -354.347
## 5000 meter buffer 3.2560 3.4670 2.7820 0.4121 6.18200 15 -355.372
## 4750 meter buffer 4.0390 2.2960 1.8190 -0.7397 4.20400 15 -355.771
## 4000 meter buffer 2.2920 2.3790 1.8270 0.1488 5.63400 15 -358.691
## 4250 meter buffer 3.2930 2.2050 1.3790 -0.4315 4.18700 15 -358.827
## 2750 meter buffer 1.4730 1.0590 1.0570 -0.1733 6.52300 15 -360.886
## 3000 meter buffer 1.3180 2.3720 1.5150 0.3454 6.65000 15 -360.909
## 2500 meter buffer 1.7490 1.0560 1.5560 0.2728 7.29400 15 -363.276
## 3250 meter buffer 1.6230 2.3810 1.8560 0.5015 5.15700 15 -363.668
## 2250 meter buffer 1.6820 0.4448 0.8810 -0.4253 5.53200 15 -369.598
## 2000 meter buffer 2.0240 0.4713 0.7796 -0.1708 4.52500 15 -371.687
## 1750 meter buffer 1.9840 0.1992 0.9361 -0.2809 3.89400 15 -374.290
## 1250 meter buffer 1.4540 -0.5564 0.4263 0.2381 -0.34050 15 -380.458
## 1500 meter buffer 1.7300 1.0630 1.0500 0.7926 1.03500 15 -382.732
## 500 meter buffer 0.8089 -1.8060 -1.4540 -0.7659 0.07479 15 -385.497
## 750 meter buffer 1.1670 -0.9409 -0.8997 -1.0550 0.67000 15 -386.228
## 1000 meter buffer 1.1790 -0.9114 -0.1012 -0.6557 -0.31860 15 -387.983
## AICc delta weight
## 3750 meter buffer 733.9 0.00 0.814
## 4500 meter buffer 737.0 3.18 0.166
## 3500 meter buffer 742.2 8.37 0.012
## 5000 meter buffer 744.3 10.42 0.004
## 4750 meter buffer 745.1 11.22 0.003
## 4000 meter buffer 750.9 17.06 0.000
## 4250 meter buffer 751.2 17.33 0.000
## 2750 meter buffer 755.3 21.45 0.000
## 3000 meter buffer 755.3 21.49 0.000
## 2500 meter buffer 760.1 26.23 0.000
## 3250 meter buffer 760.9 27.01 0.000
## 2250 meter buffer 772.7 38.87 0.000
## 2000 meter buffer 776.9 43.05 0.000
## 1750 meter buffer 782.1 48.25 0.000
## 1250 meter buffer 794.4 60.59 0.000
## 1500 meter buffer 799.0 65.14 0.000
## 500 meter buffer 804.5 70.67 0.000
## 750 meter buffer 806.0 72.13 0.000
## 1000 meter buffer 809.5 75.64 0.000
## Models ranked by AICc(x)
## Random terms (all models):
## cond(1 | array)For deer the top model was also the 3750 buffer
Let’s get the model summary
summary(deer_mods_no250$`3750 meter buffer`)## Family: binomial ( logit )
## Formula:
## cbind(`white-tailed_deer`, `absent_white-tailed_deer`) ~ seismic_lines +
## pipeline + borrowpits + wellsites + roads + trails + lc_class20 +
## lc_class34 + lc_class50 + lc_class110 + lc_class210 + lc_class220 +
## lc_class230 + (1 | array)
## Data: .
##
## AIC BIC logLik deviance df.resid
## 730.3 775.7 -350.2 700.3 137
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## array (Intercept) 2.364 1.537
## Number of obs: 152, groups: array, 4
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 6.8598 1.8961 3.618 0.000297 ***
## seismic_lines 2.1673 1.1795 1.837 0.066142 .
## pipeline 2.5130 1.6884 1.488 0.136650
## borrowpits -16.9085 5.4171 -3.121 0.001800 **
## wellsites 6.9960 2.8838 2.426 0.015268 *
## roads 2.7931 1.3925 2.006 0.044881 *
## trails 0.1252 1.2651 0.099 0.921174
## lc_class20 -6.0933 2.0278 -3.005 0.002657 **
## lc_class34 -28.0335 14.8735 -1.885 0.059457 .
## lc_class50 -9.7312 1.7469 -5.571 2.54e-08 ***
## lc_class110 -11.3860 2.2502 -5.060 4.19e-07 ***
## lc_class210 -7.0520 1.6180 -4.359 1.31e-05 ***
## lc_class220 -14.0322 2.2148 -6.336 2.36e-10 ***
## lc_class230 -5.2171 2.2354 -2.334 0.019601 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Coefficient plots
This is also a dog shit way to do this but I need to get this done
coeffs <- confint(black_bear_mods_no250$`500 meter buffer`) %>%
as_tibble() %>%
# subset to just the HFI variables for these plots
slice_head(n = 6) %>%
# add a column where we can put the feature names
rowid_to_column() %>%
# rename the columns
rename('lower' = `2.5 %`,
'upper' = `97.5 %`,
'estimate' = Estimate,
'feature' = rowid) %>%
# rename the entries to features, need to look at the order the features are in from the model summary and ensure it matches
mutate(feature = as.factor(feature),
feature = recode(feature,
'1' = 'seismic_lines',
'2' = 'pipeline',
'3' = 'borrowpits',
'4' = 'wellsites',
'5' = 'roads',
'6' = 'trails'))
# print
coeffs## # A tibble: 6 × 4
## feature lower upper estimate
## <fct> <dbl> <dbl> <dbl>
## 1 seismic_lines -3.09 -0.534 -1.81
## 2 pipeline -1.39 0.345 -0.520
## 3 borrowpits -2.43 0.0636 -1.18
## 4 wellsites -16.5 -1.27 -8.88
## 5 roads -1.20 2.18 0.492
## 6 trails -1.54 0.788 -0.373Try ad do this for all top models at once with purrr
top_models <- list(black_bear_mods_no250$`500 meter buffer`,
caribou_mods_no250$`1250 meter buffer`,
coyote_mods_no250$`4500 meter buffer`,
fisher_mods_no250$`2000 meter buffer`,
wolf_mods_no250$`4500 meter buffer`,
lynx_mods_no250$`1000 meter buffer`,
moose_mods_no250$`3750 meter buffer`,
fox_mods_no250$`3750 meter buffer`,
deer_mods_no250$`3750 meter buffer`) %>%
# use purrr to create coefficient table for all models
purrr::map(
~.x %>%
confint() %>%
as_tibble() %>%
# subset to just the HFI variables for these plots
slice_head(n = 6) %>%
# add a column where we can put the feature names
rowid_to_column() %>%
# rename the columns
rename('lower' = `2.5 %`,
'upper' = `97.5 %`,
'estimate' = Estimate,
'feature' = rowid) %>%
# rename the entries to features, need to look at the order the features are in from the model summary and ensure it matches
mutate(feature = as.factor(feature),
feature = recode(feature,
'1' = 'seismic_lines',
'2' = 'pipeline',
'3' = 'borrowpits',
'4' = 'wellsites',
'5' = 'roads',
'6' = 'trails'))) %>%
purrr::set_names('Black bear',
'Caribou',
'Coyote',
'Fisher',
'Grey wolf',
'Lynx',
'Moose',
'Red fox',
'White-tailed deer')Merge data into one data frame
coeffs_df_all <- list_rbind(top_models,
names_to = 'species') %>%
mutate(species = as.factor(species),
# add phylopic uuid for each species for plotting
# the uuid is extracted using getuuid with the species name as name = ''
uuid = case_when(species == 'Black bear' ~ get_uuid(name = 'Ursus americanus'),
species == 'Caribou' ~ get_uuid(name = 'Rangifer tarandus'),
species == 'Coyote' ~ get_uuid(name = 'Canis latrans'),
species == 'Fisher' ~ get_uuid(name = 'Pekania pennanti'),
species == 'Grey wolf' ~ get_uuid(name = 'Canis lupus'),
species == 'Lynx' ~ get_uuid(name = 'Lynx lynx'),
species == 'Moose' ~ '74eab34a-498c-4614-aece-f02361874f79',
species == 'Red fox' ~ get_uuid(name = 'Vulpes vulpes'),
species == 'White-taield deer' ~ '56f6fdb2-15d0-43b5-b13f-714f2cb0f5d0')) %>%
# need to remove problematic estimate which is going to skew plot since its so large compared to others
filter(!species == 'Red fox' &
!feature == 'wellsites')After plotting the moose image I don’t like it, let’s manually replace it in the data
# I went on the phylopic website and saw there are three images for moose, I like the last one better so we will use it
get_uuid(name = 'Alces alces',
n = 3)## [1] "df2d0ad0-adb0-49d7-afe5-edc6cad21064"
## [2] "1a20a65d-1342-4833-a9dd-1611b9fb383c"
## [3] "74eab34a-498c-4614-aece-f02361874f79"get_uuid(name = 'Odocoileus virginianus',
n = 3)## [1] "4584be20-4514-4673-a3e8-97e2a6a10e57"
## [2] "49a5a5db-047e-4e17-849b-9f96a93f0d2b"
## [3] "56f6fdb2-15d0-43b5-b13f-714f2cb0f5d0"# Then I manually copied this uuid and replaces it in the code aboveTry plotting all
ggplot(coeffs_df_all, aes(x = feature,
y = estimate,
group = uuid)) +
geom_errorbar(aes(ymin = lower,
ymax = upper,
color = feature),
width = 0.4,
linewidth = 1,
position = position_dodge(width = 1)) +
# add points for each estimate for each covariate and use position = position_dodge to shift the points so all the species don't plot on top of one another
geom_phylopic(aes(x = feature,
y = estimate,
uuid = uuid),
position = position_dodge(width = 1),
size = 4)## Warning: Removed 5 rows containing missing values.
ggplot(coeffs, aes(x = feature, y = estimate)) +
geom_point(size = 3, position = position_dodge(width = 0.5)) +
geom_errorbar(aes(ymin = lower, ymax = upper),
width = 0.4,
linewidth = 1,
position = position_dodge(width = 0.5)) +
geom_hline(yintercept = 0, linetype = "dashed")+
scale_color_manual(values = c("#56B4E9", "#009E73"), name = "Spatial Scale")+
theme_classic()+
ggtitle("Moose Response to Anthropogenic Disturbance Features")+
ylab("Coefficient Estimate \n \u00B1 95% CI")+
scale_x_discrete(labels =c("Borrowpits", "Harvest\nAreas", "Industrial\nSites", "Pipelines","Roads", "Seismic\nLines", "Trails", "Transmission\nLines"))+
theme(axis.title.x = element_blank(),
axis.text.x = element_text(size = 12),
axis.title.y = element_text(size = 14),
legend.title = element_text(size = 12),
plot.title = element_text(size = 15, hjust = 0.5))
If all else fails can use plot_model function